text:
# %% imports

using PyCall
gym = pyimport("gym")

include("../../Modules/Gridcells/StateSpikeGaussGrid.jl")
include("../../Modules/Networks/net_res.jl")
include("simviz.jl")

# %% Build the grids

include("../../Modules/Gridcells/HexGrid.jl")

env = gym.make("CartPole-v1")

cells_tmp_pos = HexModule[]
cells_tmp_vel = HexModule[]

maxs = env.observation_space.high
mins = env.observation_space.low

borpos = [convert(Float64,mins[1]),convert(Float64,maxs[1])]
borvpos = [-15.0,15.0]
#bortheta = [convert(Float64,mins[3]),convert(Float64,maxs[3])]
bortheta = [-pi,pi]
borvtheta = [-15.0,15.0]
dilats = [1, sqrt(2), 2, 2*sqrt(2), 4]

for neuron in 1:16 # parameter
    res = 0.07 # parameter
    dilat = rand()*(1.4-0.7)+0.7 # parameter
    theta = rand()*(2pi/6)-pi/6
    mod_pos = HexModule(res,dilat = dilat, theta = theta, bor_x = borpos, bor_y = bortheta)
    mod_vel = HexModule(res,dilat = dilat, theta = theta, bor_x = borvpos, bor_y = borvtheta)
    push!(cells_tmp_pos,mod_pos)
    push!(cells_tmp_vel,mod_vel)
end

const cells_pos = cells_tmp_pos
const cells_vel = cells_tmp_vel

env.close()

# %% import the trained cartpole

torch = pyimport("torch")
nn    = pyimport("torch.nn")
f     = pyimport("torch.nn.functional")

@pydef mutable struct DQLnet <: nn.Module
    function __init__(self)
        pybuiltin(:super)(DQLnet, self).__init__()
        self.fc1 = nn.Linear(4,16)
        self.fc2 = nn.Linear(16,16)
        self.fc3 = nn.Linear(16,2)
    end

    function forward(self, s)
        s = f.relu(self.fc1(s))
        s = f.relu(self.fc2(s))
        s = self.fc3(s)
        return s
    end
end

Qnet = DQLnet()
Qnet.load_state_dict(torch.load("/Users/titou/Documents/PFE/RL-STDP/Julia/CartPole/saved_cartpole.net"))

# %% Function

function count_spikes(arch::Array{Int64}, spiked::Array{Int64}, layer::Int64)
    count_layer = zeros(arch[layer])
    interval = [sum(arch[1:layer-1])+1, sum(arch[1:layer])]
    spikes_in = filter(x-> interval[1] <= x <= interval[2], spiked)
    spikes_in .-= sum(arch[1:layer-1])
    count_layer[spikes_in] .+= 1
    return count_layer
end

function left_right(count_spikes::Array{Float64}, n_left::Int64, n_right::Int64)
    sl = sum(count_spikes[1:n_left])
    sr = sum(count_spikes[n_left+1:end])
    (sl >= sr) ? action = 0 : action = 1
    return action, sl, sr
end

# %% Simulation

global s_hist = [[],[]]

function run_sim(net::Network, n_steps::Int64, render::Bool=false;
                 method::String = "snn")
    env = gym.make("CartPole-v1")
    env.seed(0)
    obs = env.reset()
    win_size = 10 # ms
    reward_count = 0
    spike_history = Array{Int64}[]
    @inbounds for j in 1:n_steps
        count_tot = zeros(net.arch[end])
        spikes = latency_spikes(obs, cells_pos, cells_vel, win_size)
        influence = 1
        @inbounds for msec in 1:win_size
            (msec==10) ? train=true : train=false
            spiked = step!(net, spikes[msec], train)
            if msec >= length(net.arch)-1
                count_tot .+= influence .* count_spikes(net.arch, spiked, length(net.arch))
                influence *= 0.95
                push!(spike_history,spiked)
            end
        end
        if method == "snn"
            n_left = div(net.arch[end],2)
            n_right = net.arch[end] - n_left
            action, sl, sr = left_right(count_tot, n_left, n_right)
            push!(s_hist[1], sl)
            push!(s_hist[2], sr)
        elseif method == "trained"
            n_left = div(net.arch[end],2)
            n_right = net.arch[end] - n_left
            action = torch.argmax(Qnet(torch.tensor(obs).float())).item()
            _, sl, sr = left_right(count_tot, n_left, n_right)
            push!(s_hist[1], sl)
            push!(s_hist[2], sr)
        elseif method == "random"
            n_left = div(net.arch[end],2)
            n_right = net.arch[end] - n_left
            action = rand(0:1)
            _, sl, sr = left_right(count_tot, n_left, n_right)
            push!(s_hist[1], sl)
            push!(s_hist[2], sr)
        end
        obs_new, reward, done, _ = env.step(action)
        reward_count += reward
        if reward == 1
            net.da += 0.1 # parameter
        elseif reward == 0
            net.da -= 0.2
        end
        if render==true
            env.render()
        end
        obs = obs_new
    end
    return spike_history, reward_count
end

# %% Create the network

# randomize the syn_weight mean = 2, std = 2 # parameters
function randomize!(net::Network)
    m = net.param["wmax"]/2
    std = net.param["wmax"]/10
    net.s_ee .= randn(net.n_exc, net.n_exc)*std .+ m
    net.s_ee .= clamp.(net.s_ee, 0.0, net.param["wmax"])
    net.s_ee .*= net.connection_ee
end

function suma(arch::Array{Int64}, i::Int64)
    return sum(arch[1:i])
end

function randomize_layer!(net::Network)
    a = net.arch
    m_base = net.param["wmax"]/2
    std_base = net.param["wmax"]/10
    con_base = a[1]*a[2]
    for l in 1:length(a)-1
        n_con = a[l]*a[l+1]
        ratio = n_con/con_base
        m = m_base/(ratio)^0.5
        std = std_base/(ratio)^0.5
        net.s_ee[suma(a,l-1)+1:suma(a,l),suma(a,l)+1:suma(a,l+1)] .= randn(a[l], a[l+1])*std .+ m
    end
    net.s_ee .= clamp.(net.s_ee, 0.0, net.param["wmax"])
    net.s_ee .*= net.connection_ee
end

parameters = YAML.load(open("/Users/titou/Documents/PFE/RL-STDP/Julia/Modules/cfg.yml"))

n = Network([128,64,64,80], parameters)
randomize_layer!(n)

global rewards = Int64[]

# %% run the simulation

@time spike_history, reward = run_sim(n, 200, method = "snn")
push!(rewards, reward)

for _ in 1:99
    @time spike_history, reward = run_sim(n, 200, method = "snn")
     push!(rewards, reward)
end

# %% Plots difference
x = collect(1:length(s_hist[1]))
y = s_hist[2] .- s_hist[1]
plot(x,y,
     xlabel = "Time step cartpole",
     ylabel = "Sum spikes left minus right",
     title = "Difference between left or right",
     legend = false)

# %% Plot s_hist
x = collect(1:length(s_hist[1]))
plot(x,s_hist,
     xlabel = "Time step cartpole",
     ylabel = "Sum spikes left and right",
     lab = ["left" "right"],
     title = "Left or Right plot")

# %% spike history plot per layer
layer_spike_hist(n.arch, spike_history, 4)

# %% spike distribution for all layers
spike_histogram(n.arch, spike_history)

# %% weights visualization heatmap
weight_viz(n, "heatmap")

# %% weights visualization histogram
weight_viz(n, "histogram")

# %% visualize network connection
 heatmap(n.connection)

# %% plot rewards
plot(rewards)

max_stars_repo_name: SuReLI/RL-STDP

max_issues_repo_name: SuReLI/RL-STDP

max_forks_repo_name: SuReLI/RL-STDP


text:
module json

end

ext: jl

max_stars_repo_name: tmaone/xsh

max_issues_repo_name: turing-complete/xsh

max_forks_repo_name: turing-complete/xsh


text:
function run_program!(p, input = Channel{Int}(), output=Channel{Int}())
    index = 0
    base = Ref(0)
    modes = zeros(Int, 3)
    program = Dict((i-1)=>v for (i,v) in enumerate(p))
    while get(program, index, 0) != 99
        index = process_opcode!(index, program, input, output, base, modes)
    end
    program
end

function process_opcode!(index, p, input, output, base, modes)
    @inbounds begin
        opcode = p[index] % 100
        modes[1], modes[2], modes[3] = p[index] ÷ 100 % 10, p[index] ÷ 1000 % 10, p[index] ÷ 10000 % 10
        if opcode == 1
            ind = modes[3] == 0 ? p[index + 3] : p[index + 3] + base[]
            p[ind] = get(p, index + 1, modes, 1, base) + get(p, index + 2, modes, 2, base)
            return index + 4
        elseif opcode == 2
            ind = modes[3] == 0 ? p[index + 3] : p[index + 3] + base[]
            p[ind] = get(p, index + 1, modes, 1, base) * get(p, index + 2, modes, 2, base)
            return index + 4
        elseif opcode == 3
            ind = modes[1] == 0 ? p[index + 1] : p[index + 1] + base[]
            p[ind] = take!(input)
            return index + 2
        elseif opcode == 4
            put!(output, get(p, index + 1, modes, 1, base))
            return index + 2
        elseif opcode == 5
            get(p, index + 1, modes, 1, base) != 0 && return get(p, index + 2, modes, 2, base)
            return index + 3
        elseif opcode == 6
            get(p, index + 1, modes, 1, base) == 0 && return get(p, index + 2, modes, 2, base)
            return index + 3
        elseif opcode == 7
            ind = modes[3] == 0 ? p[index + 3] : p[index + 3] + base[]
            p[ind] = get(p, index + 1, modes, 1, base) < get(p, index + 2, modes, 2, base) ? 1 : 0
            return index + 4
        elseif opcode == 8
            ind = modes[3] == 0 ? p[index + 3] : p[index + 3] + base[]
            p[ind] = get(p, index + 1, modes, 1, base) == get(p, index + 2, modes, 2, base) ? 1 : 0
            return index + 4
        elseif opcode == 9
            base[] += get(p, index + 1, modes, 1, base)
            return index + 2
        else
            @error("invalid opcode at index $index : $(p[index])")
            return 0
        end
    end
end

function Base.get(p::Dict{Int, Int}, i::Int, modes::Vector{Int}, ind_mode::Int, base::Ref{Int})
    m = modes[ind_mode]
    if m == 0
        # p[p[i]]
        return get(p, get(p, i, 0), 0)
    elseif m == 1
        # p[i]
        return get(p, i, 0)
    elseif m == 2
        # p[p[i] + base]
        return get(p, get(p, i, 0) + base[], 0)
    end
end


text:
abstract type Wavelet end

function minphase(wav::Array{T}) where T<:AbstractFloat
    n = length(wav)

    # zero pad
    wav = [wav; zeros(T,n)]

    # fft
    WAV = fft(wav)

    # log(|WAV|)
    logabswav = log.(abs.(WAV) .+ Base.eps(T))

    # Hilbert transform of log(|WAV|)
    hlogabswav = imag.(hilbert(logabswav)) # DSP.hilbert gives the analytic signal

    # min-phase:
    PHI = -hlogabswav

    # min-phase wavelet:
    WAV = abs.(WAV) .* exp.(im * PHI)
    wav = real.(ifft(WAV))
    wav[1:n]
end

function phaseshift(wav::Array{T}, phi::T) where T<:AbstractFloat
    if isapprox(phi, 0.0)
        return wav
    end
    n = length(wav)

    # zero pad
    wav = [wav; zeros(n)]

    # fft
    WAV = rfft(wav)

    # phase shift
    for i = 1:length(WAV)
        WAV[i] *= exp.(-im*phi)
    end
    irfft(WAV, length(wav))[1:n]
end

function costap(w::Wavelet, t::Array{T}) where T<:AbstractFloat
    if w.tmax < 0
        return ones(T,length(t))
    end
    itstrt = clamp(round(Int, .9*w.tmax/t[end]),1,length(t))
    itends = clamp(round(Int, (w.tmax/t[end])*length(t)),1,length(t))
    tap = ones(T,length(t))
    for i = 1:length(t)
        if itstrt <= i <= itends
            tap[i] = .5*(1.0 + cos(pi*(i-itstrt)/(itends-itstrt)))
        elseif i > itends
            tap[i] = 0.0
        end
    end
    tap
end

mutable struct WaveletSine <: Wavelet
    a::Float64
    f::Float64
    integrate::Bool
end
WaveletSine(;a = 1.0, f = 10.0, integrate = false) = WaveletSine(a, f, integrate)
copy(w::WaveletSine) = WaveletSine(w.a,w.f,w.integrate)

mutable struct WaveletRicker <: Wavelet
    a::Float64
    f::Float64
    tmax::Float64
    phase::Float64
    integrate::Bool
end
WaveletRicker(;a = 1.0, f = 15.0, tmax = -1.0, phase = 0.0, integrate = false) = WaveletRicker(a, f, tmax, phase, integrate)
copy(w::WaveletRicker) = WaveletRicker(w.a,w.f,w.tmax,w.phase,w.integrate)

mutable struct WaveletMinPhaseRicker <: Wavelet
    a::Float64
    f::Float64
    tmax::Float64
    integrate::Bool
end
WaveletMinPhaseRicker(;a = 1.0, f = 15.0, tmax = -1.0, integrate = false) = WaveletMinPhaseRicker(a, f, tmax, integrate)
copy(w::WaveletMinPhaseRicker) = WaveletMinPhaseRicker(w.a,w.f,w.tmax,w.integrate)

mutable struct WaveletDerivRicker <: Wavelet
    a::Float64
    f::Float64
    tmax::Float64
end
WaveletDerivRicker(;a = 1.0, f = 15.0, tmax = -1.0) = WaveletDerivRicker(a, f, tmax)
copy(w::WaveletDerivRicker) = WaveletDerivRicker(w.a,w.f,w.tmax)

mutable struct WaveletCausalRicker <: Wavelet
    a::Float64
    f::Float64
    tol::Float64
    tmax::Float64
    integrate::Bool
end
WaveletCausalRicker(;a = 1.0, f = 15.0, tol = 1.e-4, tmax = -1.0, integrate = false) = WaveletCausalRicker(a, f, tol, tmax, integrate)
copy(w::WaveletCausalRicker) = WaveletCausalRicker(w.a,w.f,w.tol,w.tmax,w.integrate)

mutable struct WaveletOrmsby <: Wavelet
    a::Float64
    f1::Float64
    f2::Float64
    f3::Float64
    f4::Float64
    tmax::Float64
    integrate::Bool
end
WaveletOrmsby(;a = 1.0, f1 = 5.0, f2 = 10.0, f3 = 40.0, f4 = 45.0, tmax = -1.0, integrate = false) = WaveletOrmsby(a, f1, f2, f3, f4, tmax, integrate)
copy(w::WaveletOrmsby) = WaveletOrmsby(w.a,w.f1,w.f2,w.f3,w.f4,w.tmax,w.integrate)

mutable struct WaveletMinPhaseOrmsby <: Wavelet
    a::Float64
    f1::Float64
    f2::Float64
    f3::Float64
    f4::Float64
    tmax::Float64
    integrate::Bool
end
WaveletMinPhaseOrmsby(;a = 1.0, f1 = 5.0, f2 = 10.0, f3 = 40.0, f4 = 45.0, tmax = -1.0, integrate = false) = WaveletMinPhaseOrmsby(a, f1, f2, f3, f4, tmax, integrate)
copy(w::WaveletMinPhaseOrmsby) = WaveletMinPhaseOrmsby(w.a,w.f1,w.f2,w.f3,w.f4,w.tmax,w.integrate)

function get(w::WaveletSine, t::Array{T}) where T<:AbstractFloat
    wav = T(w.a) .* sin.(T(2 * pi * w.f) .* t)
    w.integrate ? cumsum(wav, dims=1) : wav
end

function get(w::WaveletRicker, t::Array{T}) where T<:AbstractFloat
    wav = convert(Array{T},w.a * (1.0 .- 2.0 * (pi * w.f * t).^2) .* exp.(-(pi * w.f * t).^2))
    wav = phaseshift(wav, T(w.phase))
    wav .*= costap(w,t)
    w.integrate ? cumsum(wav, dims=1) : wav
end

function get(w::WaveletMinPhaseRicker, t::Array{T}) where T<:AbstractFloat
    wav = minphase(get(WaveletRicker(w.a,w.f,-1.0,0.0,false),t .- mean(t)))
    wav .*= costap(w,t)
    w.integrate ? cumsum(wav, dims=1) : wav
end

function get(w::WaveletDerivRicker, t::Array{T}) where T<:AbstractFloat
    wav1 = get(WaveletRicker(w.a,w.f,-1.0,0.0,false),t)
    wav = similar(wav1)
    wav[1] = wav1[1]
    for i = 2:length(wav1)
        wav[i] = wav1[i] - wav1[i-1]
    end
    wav .* costap(w,t)
end

function get(w::WaveletCausalRicker, t::Array{T}) where T<:AbstractFloat
    wav = get(WaveletRicker(w.a,w.f,-1.0,0.0,false),t .- mean(t))
    idx = findall(ww->abs(ww)>w.tol*w.a,wav)[1]
    rng = idx:length(wav)
    wav[1:length(rng)] = wav[rng]
    wav[length(rng)+1:end] .= 0.0
    wav .*= costap(w,t)
    w.integrate ? cumsum(wav, dims=1) : wav
end

function get(w::WaveletOrmsby, t::Array{T}) where T<:AbstractFloat
    numerator(f,t) = sinc.(f*t).^2 .* (pi*f)^2
    pf43 = T(pi*(w.f4-w.f3))
    pf21 = T(pi*(w.f2-w.f1))
    wav = numerator(T(w.f4),t)/pf43 - numerator(T(w.f3),t)/pf43 - numerator(T(w.f2),t)/pf21 + numerator(T(w.f1),t)/pf21
    wav ./= maximum(abs, wav)
    wav .*= w.a
    wav .-= mean(wav)
    wav = wav .* costap(w,t)
    w.integrate ? cumsum(wav, dims=1) : wav
end

function get(w::WaveletMinPhaseOrmsby, t::Array{T}) where T<:AbstractFloat
    wav = minphase(get(WaveletOrmsby(w.a,w.f1,w.f2,w.f3,w.f4,w.tmax,false),t .- mean(t)))
    wav ./= maximum(abs, wav)
    wav .*= w.a
    wav .-= mean(wav)
    wav .*= costap(w,t)
    w.integrate ? cumsum(wav, dims=1) : wav
end

max_stars_repo_name: ChevronETC/WaveFD

max_issues_repo_name: ChevronETC/WaveFD

max_forks_repo_name: ChevronETC/WaveFD


text:
using ModelingToolkit
using Test

@parameters t
@variables x(t) y(t) z(t)

null_op = 0*t
@test isequal(simplify_constants(null_op), 0)

one_op = 1*t
@test isequal(simplify_constants(one_op), t)

identity_op = Operation(identity,[x])
@test isequal(simplify_constants(identity_op), x)

minus_op = -x
@test isequal(simplify_constants(minus_op), -1*x)
simplify_constants(minus_op)


text:
#=

# Package options

Certain features of this library are configurable via [Preferences.jl](https://github.com/JuliaPackaging/Preferences.jl).

|    Preference     |              Description              |  Default  |
|:-----------------:|:-------------------------------------:|:---------:|
| `LOG_DESTRUCTION` | Log the destruction of Vulkan handles | `"false"` |
| `USE_DISPATCH_TABLE` | Retrieve and store function pointers in a [dispatch table](@ref Dispatch) to speed up API calls and facilitate the use of extensions | `"true"` |

## Destruction logging

To debug the reference counting mechanism used through finalizers for handle destruction, it is possible to print when a finalizer is run and the resulting effect (freeing the object, or simply decrementing a reference counter). To enable this, set the preference `LOG_DESTRUCTION` to `"true"`.

=#

max_stars_repo_name: JuliaGPU/VulkanAbstraction

max_issues_repo_name: JuliaGPU/VulkanAbstraction

max_forks_repo_name: JuliaGPU/VulkanAbstraction


text:
cd("C:/Users/Philippe/Documents/Research/Irregular Coulomb Potential")
using SincFun, DESincEig, PyPlot, Distributions
include("C:/Users/Philippe/Documents/Julia/Iregular Potential/Iregular_V1.jl")
Range=collect(1:1:40)
n = 1
@time (RESULTS, All_Abs_Error_Approx , hoptimal , n, MatrixSizes, All_Eigenvectors, MN,ap,an) = Laurent_cases(1,Range=Range,tol=5e-10,enum=(n,NaN))

savefig("Potential3_1.eps")


text:
"""
GenX: An Configurable Capacity Expansion Model
Copyright (C) 2021,  Massachusetts Institute of Technology
This program is free software; you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation; either version 2 of the License, or
(at your option) any later version.
This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
GNU General Public License for more details.
A complete copy of the GNU General Public License v2 (GPLv2) is available
in LICENSE.txt.  Users uncompressing this from an archive may not have
received this license file.  If not, see <http://www.gnu.org/licenses/>.
"""

@doc raw"""
	write_storagedual(path::AbstractString, inputs::Dict, setup::Dict, EP::Model)

Function for reporting dual of storage level (state of charge) balance of each resource in each time step.
"""
function write_storagedual(path::AbstractString, inputs::Dict, setup::Dict, EP::Model)
	dfGen = inputs["dfGen"]

	G = inputs["G"]     # Number of resources (generators, storage, DR, and DERs)
	T = inputs["T"]     # Number of time steps (hours)
	Z = inputs["Z"]     # Number of zones

	START_SUBPERIODS = inputs["START_SUBPERIODS"]
	INTERIOR_SUBPERIODS = inputs["INTERIOR_SUBPERIODS"]
	REP_PERIOD = inputs["REP_PERIOD"]

	hours_per_subperiod = inputs["hours_per_subperiod"] #total number of hours per subperiod

	# # Dual of storage level (state of charge) balance of each resource in each time step
	dfStorageDual = DataFrame(Resource = inputs["RESOURCES"], Zone = dfGen[!,:Zone])
	# Define an empty array
	x1 = Array{Float64}(undef, G, T)
	dual_values =Array{Float64}(undef, G, T)

	# Loop over W separately hours_per_subperiod
	for y in 1:G
		if y in inputs["STOR_ALL"]
			print(y)
			if setup["OperationWrapping"]==1 && y in inputs["STOR_LONG_DURATION"]
				for w in 1:REP_PERIOD
					x1[y,hours_per_subperiod*(w-1)+1] = dual.(EP[:cSoCBalLongDurationStorageStart][w,y])
				end
			else
				for t in START_SUBPERIODS
					x1[y,t] = dual.(EP[:cSoCBalStart][t,y])
				end
			end
			for t in INTERIOR_SUBPERIODS
				x1[y,t] = dual.(EP[:cSoCBalInterior][t,y]) #Use this for getting dual values and put in the extracted codes from PJM
			end
		else
			x1[y,:] = zeros(T,1) # Empty values for the resource with no ability to store energy
		end
	end

	# Incorporating effect of time step weights (When OperationWrapping=1) and Parameter scaling (ParameterScale=1) on dual variables
	for y in 1:G
		if setup["ParameterScale"]==1
			dual_values[y,:] = x1[y,:]./inputs["omega"] *ModelScalingFactor
		else
			dual_values[y,:] = x1[y,:]./inputs["omega"]
		end
	end

	dfStorageDual=hcat(dfStorageDual, DataFrame(dual_values, :auto))
	rename!(dfStorageDual,[Symbol("Resource");Symbol("Zone");[Symbol("t$t") for t in 1:T]])

	CSV.write(joinpath(path, "storagebal_duals.csv"), dftranspose(dfStorageDual, false), writeheader=false)
end

max_stars_repo_name: Betristor/LTESOM

max_issues_repo_name: Betristor/LTESOM

max_forks_repo_name: Betristor/LTESOM


text:
cellmerge(x::Number, y::Number) = x ≈ y ? x : contradiction

using IntervalArithmetic

function cellmerge(x::Interval, y::Interval)
  xy = intersect(x, y)
  isempty(xy) ? contradiction : xy
end


text:
### A Pluto.jl notebook ###
# v0.12.20

using Markdown
using InteractiveUtils

# This Pluto notebook uses @bind for interactivity. When running this notebook outside of Pluto, the following 'mock version' of @bind gives bound variables a default value (instead of an error).
macro bind(def, element)
    quote
        local el = $(esc(element))
        global $(esc(def)) = Core.applicable(Base.get, el) ? Base.get(el) : missing
        el
    end
end

# ╔═╡ 747d446a-dfeb-11ea-3533-c9404fd41688
begin
	using IndividualDisplacements, MeshArrays, DataFrames
	using Statistics, PlutoUI, StatsPlots, OrdinaryDiffEq

	p=dirname(pathof(IndividualDisplacements))
    include(joinpath(p,"../examples/helper_functions.jl"))
    #include(joinpath(p,"../examples/recipes_plots.jl"))
	𝑃,𝐷,Γ=OCCA_FlowFields()
	tmp=(Γ = Γ, m = "OCCA")
    𝐷=merge(𝐷,tmp)
	✓ = "😃"
	"$✓ Set up packages, gridded domain, etc"
end

# ╔═╡ bf19d29c-e70e-11ea-0153-d3a49981d56c
md"""longitude,latitude of South West corner =
$(@bind lon0 NumberField(-180.0:10:180.0; default=0.0))
,
$(@bind lat0 NumberField(-90.0:5:90.0; default=-50.0))
"""

# ╔═╡ 4935fd46-e70f-11ea-386c-f9c444a20644
md"""depth level (integer), set size (nb of particles) =
$(@bind klev NumberField(1:10; default=5))
,
$(@bind npar NumberField(50:50:1000; default=100))
"""

# ╔═╡ 9c80e722-e70f-11ea-22a6-0be2e85f3b8b
md"""time step (for plotting)=
$(@bind tt Slider(1:26; default=13, show_value=true))
"""

# ╔═╡ e25eee9e-dfee-11ea-2a4c-3946ccb63876
begin
	lo0=lon0; lo1=lo0+5
	la0=lat0; la1=la0+5
	z_init=klev
	n_part=npar
	(particles=n_part,longitudes=(lo0,lo1),latitudes=(la0,la1),level=z_init,plot_time=tt)
end

# ╔═╡ f75fae30-dfee-11ea-18ef-259321acfa2f
begin
	lon=lo0 .+(lo1-lo0).*rand(n_part)
	lat=la0 .+(la1-la0).*rand(n_part)
	(_,_,_,_,f,x,y)=InterpolationFactors(𝐷.Γ,lon,lat)
    m=findall( (f.!==0).*((!isnan).(x)) )
	df=DataFrame(x=x[m],y=y[m],f=f[m])

	custom🔴 = DataFrame(ID=Int[], fid=Int[], x=Float32[], y=Float32[],
		z=Float32[], year=Float32[], t=Float32[])
	function custom🔧(sol,𝐹::𝐹_MeshArray3D;id=missing,𝑇=missing)
		df=postprocess_MeshArray(sol,𝐹,id=id,𝑇=𝑇)
		z=[sol[1,i,j][3] for i in 1:size(sol,2), j in 1:size(sol,3)]
		df.z=z[:]
		df.year=df.t ./86400/365
		return df
	end
	custom∫(prob) = solve(prob,Tsit5(),reltol=1e-8,abstol=1e-8,saveat=365/12*86400.0)

	𝐼=Individuals(𝑃,df.x,df.y,fill(z_init,n_part),df.f,
		(🔴=custom🔴,🔧=custom🔧,∫=custom∫))

	"$✓ Set up Individuals"
end

# ╔═╡ 9ffe84c0-dff0-11ea-2726-8924892df73a
begin
	𝑇=(0.0,𝐼.𝑃.𝑇[2])
	∫!(𝐼,𝑇)

	add_lonlat!(𝐼.🔴,𝐷.XC,𝐷.YC)

	🔴_by_ID = groupby(𝐼.🔴, :ID)
	🔴_by_t = groupby(𝐼.🔴, :t)

	ff(x) = x[1]
	tmp1=combine(🔴_by_ID,:lon => ff => :lo,:lat => ff => :la)
	𝐼.🔴.dlon=𝐼.🔴.lon - tmp1.lo[𝐼.🔴.ID]
	𝐼.🔴.dlat=𝐼.🔴.lat - tmp1.la[𝐼.🔴.ID]

	nt=length(unique(𝐼.🔴.t))

	"$✓ ∫!(𝐼,𝑇) etc"
end

# ╔═╡ f65ddffa-e63a-11ea-34a6-2fa9284e98fa
begin
	mx=20.0

	#f(x)=x[tt]-x[1]
	#f(x)=last(x).-first(x)
	#cdf = combine(🔴_by_ID,nrow,:lat => f => :dlat,:lon => f => :dlon)

	g(x)=x[tt]
	cdf = combine(🔴_by_ID,:dlat => g => :dlat,:dlon => g => :dlon)

	plt_hist=histogram2d(cdf.dlon,cdf.dlat,nbins = (10, 10),colorbar=false)
	#scatter(cdf.dlon,cdf.dlat,xlims=(-mx,mx),ylims=(-mx,mx))
	"$✓ dlon, dlat histogram2d"
end

# ╔═╡ 7d52252e-e006-11ea-2632-df2af831b52f
begin
	xx=vec(𝐷.Γ["XC"][1][:,1])
	yy=vec(𝐷.Γ["YC"][1][1,:])
	zz=transpose(log10.(𝐷.Γ["Depth"][1]))
	𝐵=(x = xx, y = yy, z = zz)

	𝐶(g::ColorGradient) = RGB[g[z] for z=LinRange(0,1,length(🔴_by_t))]
	𝐶(t::Int) = 𝐶(cgrad(:inferno))[t]
	"$✓ Set up plotting"
end

# ╔═╡ 0b12cf52-e6e3-11ea-1a01-dd0c49c9e641
begin
	plt_dlat = @df 🔴_by_t[1] density(:dlat, leg = :none, colour = 𝐶(1), ylims=(0,0.5))
	[@df 🔴_by_t[tt] density!(plt_dlat,:dlat, leg = :none, colour = 𝐶(tt)) for tt in 2:length(🔴_by_t)];
	density!(plt_dlat,🔴_by_t[tt].dlat, leg = :none, colour = :cyan, linewidth=4)
	"$✓ dlat density"
end

# ╔═╡ 6f70033a-e6cc-11ea-373e-6dcbaaa53d15
begin
	#clims=extrema(:t)
	#@df 𝐼.🔴 density(:lon, group = (:t), leg = :none, palette = cgrad(:ice))
	plt_dlon = @df 🔴_by_t[1] density(:dlon, leg = :none, colour = 𝐶(1), ylims=(0,0.5))
	[@df 🔴_by_t[tt] density!(plt_dlon,:dlon, leg = :none, colour = 𝐶(tt)) for tt in 2:length(🔴_by_t)];
	density!(plt_dlon,🔴_by_t[tt].dlon, leg = :none, colour = :cyan, linewidth=4)
	"$✓ dlon density"
end

# ╔═╡ a13d6ea6-dff1-11ea-0713-cb235e28cf79
begin
	plt_map=contourf(xx,yy,zz,clims=(-.5,4.),c = :ice,
		colorbar=false, xlims=(-180.0,180.0),ylims=(-90.0,90.0))
	scatter!(plt_map,🔴_by_t[1].lon,🔴_by_t[1].lat,c=:gold,leg=:none,
		markersize=2.0, marker = (:dot, stroke(0)) )
	scatter!(plt_map,🔴_by_t[tt].lon,🔴_by_t[tt].lat,c=:red,leg=:none,
		markersize=2.0, marker = (:dot, stroke(0)) )

	plot(plt_map,plt_hist,plt_dlon,plt_dlat)
end

# ╔═╡ Cell order:
# ╟─bf19d29c-e70e-11ea-0153-d3a49981d56c
# ╟─4935fd46-e70f-11ea-386c-f9c444a20644
# ╟─9c80e722-e70f-11ea-22a6-0be2e85f3b8b
# ╟─e25eee9e-dfee-11ea-2a4c-3946ccb63876
# ╟─a13d6ea6-dff1-11ea-0713-cb235e28cf79
# ╟─f65ddffa-e63a-11ea-34a6-2fa9284e98fa
# ╟─0b12cf52-e6e3-11ea-1a01-dd0c49c9e641
# ╟─6f70033a-e6cc-11ea-373e-6dcbaaa53d15
# ╟─9ffe84c0-dff0-11ea-2726-8924892df73a
# ╟─7d52252e-e006-11ea-2632-df2af831b52f
# ╟─f75fae30-dfee-11ea-18ef-259321acfa2f
# ╟─747d446a-dfeb-11ea-3533-c9404fd41688


text:
using BinDeps
using Compat

@BinDeps.setup

libogg = library_dependency("libogg", aliases = ["libogg"])

@static if is_apple()
  using Homebrew
  provides( Homebrew.HB, "libogg", libogg, os = :Darwin )
end

provides( AptGet, "libogg0", libogg )

@compat @BinDeps.install Dict(:libogg => :libogg)


text:
using ExpectationMaximization
using Distributions
using Test

@testset "ExpectationMaximization.jl" begin
    N = 50_000
    θ₁ = 10
    θ₂ = 5
    α = 0.2
    β = 0.3
    rtol = 5e-2
    mix_true = MixtureModel([Exponential(θ₁), Gamma(α, θ₂)], [β, 1 - β])
    y = rand(mix_true, N)
    mix_guess = MixtureModel([Exponential(1), Gamma(0.5, 1)], [0.5, 1 - 0.5])
    mix_mle = fit_mle(mix_guess, y; display = :iter, tol = 1e-3, robust = false, infos = false)

    p = params(mix_mle)[1]
    @test isapprox(β, probs(mix_mle)[1]; rtol = rtol)
    @test isapprox(θ₁, p[1]...; rtol = rtol)
    @test isapprox(α, p[2][1]; rtol = rtol)
    @test isapprox(θ₂, p[2][2]; rtol = rtol)
end


text:
module Tasks

using SearchLight, Nullables, SearchLight.Validation, TasksValidator

export Task

mutable struct Task <: AbstractModel
  ### INTERNALS
  _table_name::String
  _id::String
  _serializable::Vector{Symbol}

  ### FIELDS
  id::DbId
  content::String # 追加
  done::Bool      # 追加

  ### VALIDATION
  # validator::ModelValidator

  ### CALLBACKS
  # before_save::Function
  # after_save::Function
  # on_save::Function
  # on_find::Function
  # after_find::Function

  ### SCOPES
  # scopes::Dict{Symbol,Vector{SearchLight.SQLWhereEntity}}

  ### constructor
  Task(;
    ### FIELDS
    id = DbId(),
		content = "", # 追加
    done = false  # 追加

    ### VALIDATION
    # validator = ModelValidator([
    #   ValidationRule(:title, TasksValidator.not_empty)
    # ]),

    ### CALLBACKS
    # before_save = (m::Todo) -> begin
    #   @info "Before save"
    # end,
    # after_save = (m::Todo) -> begin
    #   @info "After save"
    # end,
    # on_save = (m::Todo, field::Symbol, value::Any) -> begin
    #   @info "On save"
    # end,
    # on_find = (m::Todo, field::Symbol, value::Any) -> begin
    #   @info "On find"
    # end,
    # after_find = (m::Todo) -> begin
    #   @info "After find"
    # end,

    ### SCOPES
    # scopes = Dict{Symbol,Vector{SearchLight.SQLWhereEntity}}()

  ) = new("tasks", "id", Symbol[],                                                ### INTERNALS
          id, content, done                                                                           ### FIELDS
          # validator,                                                                  ### VALIDATION
          # before_save, after_save, on_save, on_find, after_find                       ### CALLBACKS
          # scopes                                                                      ### SCOPES
          )
end

# 追加
function seed()
  SampleTasks = [
    ("塩を買う", true),
    ("簡易書留を出す", false),
    ("本を5冊読む", false)
  ]
  for task in SampleTasks
    Task(content = task[1], done = task[2]) |> SearchLight.save!
  end
end

end

max_stars_repo_name: yshutaro/GenieTodoList

max_issues_repo_name: yshutaro/GenieTodoList

max_forks_repo_name: yshutaro/GenieTodoList


text:
# Matrices and Vectors for the Gradients

function calctheta( z::Array{T,1}, c::SLiSeFilters.Config )
    t = zeros(Complex{Float64},length(z),1);
    for k=1:length(z)
        t[k] = theta_k( z[k], c );
    end
    return t
end

function calcthetaP( z, c::SLiSeFilters.Config )
    t = zeros(Complex{Float64},length(z),1);
    for k=1:length(z)
        t[k] = thetaP_k( z[k], c::SLiSeFilters.Config );
    end
    return t
end

# z is Hermitian
function calcZ( z, c )
    Z = zeros(Complex{Float64},length(z), length(z));
    for k=1:length(z)
        for l=1:k
            Z[k,l] = Z_kl( z[k], z[l], c );
        end
    end
    Z = Z + tril(Z, -1)';
end

# y is real symmetric
function calcY( z, c::SLiSeFilters.Config )
    Y = zeros(Complex{Float64},length(z), length(z));
    for k=1:length(z)
        for l=1:k
            Y[k,l] = Y_kl( z[k], z[l], c );
        end
    end
    Y = Y + transpose(tril(Y, -1));
end

# W is real symmetric
function  calcW( z, c::SLiSeFilters.Config )
    W = zeros(Complex{Float64},length(z), length(z));
    for k=1:length(z)
        for l=1:k
            if k == l
                W[k,l] = W_kk( z[k], c );
            else
                W[k,l] = W_kl( z[k], z[l], c );
            end
        end
    end
    W = W + transpose(tril(W, -1));
end

# W is real symmetric
function calcWtp( z, c::SLiSeFilters.Config )
    W = zeros(Complex{Float64},length(z), length(z));
    for k=1:length(z)
        for l=1:k
            if k == l
                W[k,l] = Wtp_kk( z[k], c );
            else
                W[k,l] = Wtp_kl( z[k], z[l], c );
            end
        end
    end
    W = W + transpose(tril(W, -1));
end

# X is hermitian
function calcX( z, c::SLiSeFilters.Config )
    X = zeros(Complex{Float64},length(z), length(z));
    for k=1:length(z)
        for l=1:k
            X[k,l] = X_kl( z[k], z[l], c );
        end
    end
    X = X + tril(X, -1)';
end

# Xtp is hermitian
function calcXtp( z, c::SLiSeFilters.Config )
    X = zeros(Complex{Float64},length(z), length(z));
    for k=1:length(z)
        for l=1:k
            X[k,l] = Xtp_kl( z[k], z[l], c );
        end
    end
    X = X + tril(X, -1)';
end

function theta_k( z1, c::SLiSeFilters.Config )
    case1( conj(z1), c );
end

function thetaP_k( z1, c::SLiSeFilters.Config )
    x = case1( -conj(z1), c );
end

function Wtp_kk( z1, c::SLiSeFilters.Config )
    case2( -z1, c );
end

function W_kk( z1, c::SLiSeFilters.Config )
    case2( conj(z1), c );
end

function Z_kl( z1, z2, c::SLiSeFilters.Config )
    case3( conj(z1), -z2, c );
end

function  Y_kl( z1, z2, c::SLiSeFilters.Config )
    case3( conj(z1), -conj(z2), c );
end

# 1 / (t-conj(z1)) / (t-conj(z2)) -> 1 / (t-z1) / (t-z2)
function  W_kl( z1, z2, c::SLiSeFilters.Config )
    x = case3( conj(z1), conj(z2), c );
end

# 1 / (t-conj(z1)) / (t-conj(z2)) -> 1 / (t-z1) / (t-z2)
function  Wtp_kl( z1, z2, c::SLiSeFilters.Config )
    x = case3( -z1, -z2, c );
end

#        % 1 / (t-conj(z1)) / (t-z2) -> 1 / (t-z1) / (t-z2)
function  X_kl( z1, z2, c::SLiSeFilters.Config )
    x = case3( conj(z1), z2, c );
end

#        % 1 / (t+z1)) / (t+conj(z2)) -> 1 / (t-z1) / (t-z2)
function Xtp_kl( z1, z2, c::SLiSeFilters.Config )
    x = case3( -z1, -conj(z2), c );
end

############################################################
function calcddW( z, c::SLiSeFilters.Config )
    W = zeros(Complex{Float64},length(z),length(z));
    for k = 1:length(z)
        for l = 1:length(z)
            if z[k] == z[l]
                W[k,k] = case7( z[k], c );
            else
                W[k,l] = case6( z[k], z[l], c );
            end
        end
    end
    return W
end

function calcddX( z, c::SLiSeFilters.Config )
    W = zeros(Complex{Float64},length(z),length(z));
    for k = 1:length(z)
        for l = 1:length(z)
            zk = conj(z[k]);
            zl = z[l];
            if zk == zl
                W[k,l] = case7( zk, c );
            else
                W[k,l] = case6( zk, zl, c );
            end
        end
    end
    return W
end

function calcddY( z, c::SLiSeFilters.Config )
    W = zeros(Complex{Float64},length(z),length(z));
    for k = 1:length(z)
        for l = 1:length(z)
            zk = z[k];
            zl = -z[l];
            if zk == zl
                W[k,l] = -case7( zk, c );
            else
                W[k,l] = -case6( zk, zl, c );
            end
        end
    end
    return W
end

function calcddZ( z, c::SLiSeFilters.Config )
    W = zeros(Complex{Float64},length(z),length(z));
    for k = 1:length(z)
        for l = 1:length(z)
            zk = conj(z[k]);
            zl = -z[l];
            if zk == zl
                W[k,l] = -case7( zk, c );
            else
                W[k,l] = -case6( zk, zl, c );
            end
        end
    end
    return W
end

function calcDeltaX( z, c::SLiSeFilters.Config )
    X = zeros(Complex{Float64},length(z),length(z));
    for k=1:length(z)
        for l=1:length(z)
            X[k,l] = case4( z[l], conj(z[k]), c );
        end
    end
    return X
end

function calcDeltaXp( z, c::SLiSeFilters.Config )
    X = zeros(Complex{Float64},length(z),length(z));
    for k=1:length(z)
        for l=1:length(z)
            X[k,l] = -case4( -z[l], -conj(z[k]), c );
        end
    end
    return X
end

function calcDeltaYc( z, c::SLiSeFilters.Config )
    X = zeros(Complex{Float64},length(z),length(z));
    for k=1:length(z)
        for l=1:length(z)
            X[k,l] = -case4( -z[l], z[k], c );
        end
    end
    return X
end

function calcDeltaYcp( z, c::SLiSeFilters.Config )
    X = zeros(Complex{Float64},length(z),length(z));
    for k=1:length(z)
        for l=1:length(z)
            X[k,l] = case4( z[l], -z[k], c );
        end
    end
    return X
end

function calcDeltaWc( z, c::SLiSeFilters.Config )
    X = zeros(Complex{Float64},length(z),length(z));
    for k=1:length(z)
        for l=1:length(z)
            if k == l
                X[k,l] = case5( z[l], c );
            else
                X[k,l] = case4( z[l], z[k], c );
            end
        end
    end
    return X
end

function calcDeltaWcp( z, c::SLiSeFilters.Config )
    X = zeros(Complex{Float64},length(z),length(z));
    for k=1:length(z)
        for l=1:length(z)
            if k == l
                X[k,l] = -case5( -z[l], c );
            else
                X[k,l] = -case4( -z[l], -z[k], c );
            end
        end
    end
    return X
end

function calcDeltaZ( z, c::SLiSeFilters.Config )
    X = zeros(Complex{Float64},length(z),length(z));
    for k=1:length(z)
        for l=1:length(z)
            X[k,l] = -case4( -z[l], conj(z[k]), c );
        end
    end
    return X
end

function calcDeltaZp( z, c::SLiSeFilters.Config )
    X = zeros(Complex{Float64},length(z),length(z));
    for k=1:length(z)
        for l=1:length(z)
            X[k,l] = case4( z[l], -conj(z[k]), c );
        end
    end
    return X
end

function calcDeltaThetac( z, c::SLiSeFilters.Config )
    X = zeros(Complex{Float64},length(z),1);
    for k=1:length(z)
        X[k] = case2_h( z[k], c );
    end
    return X
end

function calcDeltaTheta( z, c::SLiSeFilters.Config )
    X = zeros(Complex{Float64},length(z),1);
    for k=1:length(z)
        X[k] = case2_h( conj(z[k]), c );
    end
    return X
end

function calcDeltaThetacp( z, c::SLiSeFilters.Config )
    X = zeros(Complex{Float64},length(z),1);
    for k=1:length(z)
        X[k] = -case2_h( -z[k], c );
    end
    return X
end


text:
module CUSOLVER

using ..CuArrays
using ..CuArrays: active_context, _getindex, unsafe_free!, @argout, @workspace

using ..CUBLAS: cublasFillMode_t, cublasOperation_t, cublasSideMode_t, cublasDiagType_t
using ..CUSPARSE: cusparseMatDescr_t

using CUDAapi

using CUDAdrv
using CUDAdrv: CUstream

import CUDAnative

using CEnum

const libcusolver = Ref("libcusolver")

# core library
include("libcusolver_common.jl")
include("error.jl")
include("libcusolver.jl")

# low-level wrappers
include("util.jl")
include("wrappers.jl")

# high-level integrations
include("linalg.jl")

const _dense_handles = Dict{CuContext,cusolverDnHandle_t}()
const _dense_handle = Ref{cusolverDnHandle_t}(C_NULL)
const _sparse_handles = Dict{CuContext,cusolverSpHandle_t}()
const _sparse_handle = Ref{cusolverSpHandle_t}(C_NULL)

function dense_handle()
    if _dense_handle[] == C_NULL
        CUDAnative.maybe_initialize("CUSOLVER")
        _dense_handle[] = get!(_dense_handles, active_context[]) do
            context = active_context[]
            handle = cusolverDnCreate()
            atexit(()->CUDAdrv.isvalid(context) && cusolverDnDestroy(handle))
            handle
        end
    end
    return _dense_handle[]
end

function sparse_handle()
    if _sparse_handle[] == C_NULL
        CUDAnative.maybe_initialize("CUSOLVER")
        _sparse_handle[] = get!(_sparse_handles, active_context[]) do
            context = active_context[]
            handle = cusolverSpCreate()
            atexit(()->CUDAdrv.isvalid(context) && cusolverSpDestroy(handle))
            handle
        end
    end
    return _sparse_handle[]
end

end


text:
using TuringModels

Turing.setadbackend(:reverse_diff);
#nb Turing.turnprogress(false);

# Can't really set a Uniform[-Inf,Inf] on σ

# Turing model
@model m8_4(y) = begin
    α₁ ~ Uniform(-3000, 1000)
    α₂ ~ Uniform(-1000, 3000)
    σ ~ Truncated(Cauchy(0,1), 0, Inf)

    for i ∈ 1:length(y)
        y[i] ~ Normal(α₁ + α₂, σ)
    end
end

# Observations

y = rand(Normal(0,1), 100);

# Sample

chns = sample(m8_4(y), Turing.NUTS(0.95), 2000);

# Results rethinking

m84rethinking = "
         mean      sd     5.5%   94.5% n_eff Rhat
 a1    -861.15 558.17 -1841.89  -31.04     7 1.43
 a2     861.26 558.17    31.31 1842.00     7 1.43
 sigma    0.97   0.07     0.89    1.09     9 1.17
";

# Describe the posterior samples

describe(chns)

# End of `08/m8.4t.jl`


text:
using SafeTestsets

@time begin
    # Important: "Current Loaded Data" testset must be called first. It relies on initial null CURRENT
    @time @safetestset "Current Loaded Data" begin include("current_test.jl") end

    @time @safetestset "Save/Load Data" begin include("data_test.jl") end

    @time @safetestset "Walk Logic Tests (using Arlett)" begin include("arlett_walk_logic_test.jl") end

    @time @safetestset "Arlett Backend Tests" begin include("arlett_backend_test.jl") end
    @time @safetestset "Arlett Frontend Tests" begin include("arlett_frontend_test.jl") end

    @time @safetestset "Enz Backend Tests" begin include("enz_backend_test.jl") end
    @time @safetestset "Enz Frontend Tests" begin include("enz_frontend_test.jl") end

    @time @safetestset "PPD Backend Tests" begin include("ppd_backend_test.jl") end
    @time @safetestset "PPD Frontend Tests" begin include("ppd_frontend_test.jl") end

    @time @safetestset "Quinto Backend Tests" begin include("quinto_backend_test.jl") end
end

max_stars_repo_name: s-y-wu/HMCResearchRandomWalks

max_issues_repo_name: s-y-wu/HMCResearchRandomWalks

max_forks_repo_name: s-y-wu/HMCResearchRandomWalks


text:
"""
    function residual_demand(sys::System, res::PowerSimulations.OperationsProblemResults)

Calculates the residual demand of the slack generator for a solved OPF case. The residual
demand is defined as the remaining demand to be satisfied by the slack generator while
knowing the power output of all the other generators [1]. It can be calculated as:

    residual_slack = sum{Load[i], i in all nodes} - sum{Pg[j], j in all gens, j ≠ slack}

[1] J. Portela González, A. Muñoz San Roque, E. F. Sánchez-Úbeda, J. García-González and R.
    González Hombrados, "Residual Demand Curves for Modeling the Effect of Complex Offering
    Conditions on Day-Ahead Electricity Markets," in IEEE Transactions on Power Systems,
    vol. 32, no. 1, pp. 50-61, Jan. 2017, doi: 10.1109/TPWRS.2016.2552240.

#Arguments
- `sys::System`:                            Power system in p.u. (from PowerSystems.jl)
- `res::PowerSimulations.
    OperationsProblemResults`:              Results of the solved OPF for the system
                                             (from PowerSimulations.jl)

#Keywords
- `dual_lines_tol::Float64 = 1e-1`:         Tolerance to identify any binding lines
- `dual_gen_tol::Float64 = 1e-1`:           Tolerance to identify any binding generators

"""
function residual_demand(sys::System, res::PowerSimulations.OperationsProblemResults)

    # Define elements
    allPg_butslack = 0
    total_Pload = 0.0

    # Get all thermal generators
    gens_thermal = get_components(ThermalStandard, sys)
    all_PGenThermal = get_variables(res)[:P__ThermalStandard] #Optimised PGen

    # Get all loads of the system
    loads = get_components(PowerLoad,sys)

    # Save the active power output of all gens but the slack
    for gen_thermal in gens_thermal
        if get_bustype(gen_thermal.bus) ≠ BusTypes.REF
            allPg_butslack = allPg_butslack + all_PGenThermal[1, gen_thermal.name]
        end
    end

    # Calculate total active power load of the system
    for load in loads
        total_Pload = total_Pload + load.active_power
    end

    # Calculate the residual demand for the slack generator
    residual_slack = total_Pload - allPg_butslack

    return residual_slack
end

max_stars_repo_name: invenia/MaximumGeneratorProfitTCRDD

max_issues_repo_name: invenia/MaximumGeneratorProfitTCRDD

max_forks_repo_name: invenia/MaximumGeneratorProfitTCRDD


text:
function mese1()

    rfex = RF([π/2], [0], 0, 0, Gradient(0, 0, 0))
    rfref = RF([π], [-π/2], 0, 0, Gradient(0, 0, 0))
    mese! = MESEBlochSim(3000, 10, 32, rfex, rfref)
    s1 = Spin(1, 1000, 100, 0)
    s2 = SpinMC(1, (1, 0), (1000, 100), (100, 10), (0, 0), (Inf, Inf))
    M1 = mese!(s1)
    M2 = mese!(s2)

    show(devnull, mese!)
    show(devnull, "text/plain", mese!)
    return all(M1[i] ≈ M2[i][1] for i = 1:32)

end

@testset "MESE" begin

    @test mese1()

end


text:
function __init__()
    register(DataDep(
        "AAindex",
        """
        Dataset: AAindex
        Website: https://www.genome.jp/aaindex/

        Please cite the following references when making use of the database:

        Nakai, K., Kidera, A., and Kanehisa, M.;  Cluster analysis of
           amino acid indices for prediction of protein structure and
           function.  Protein Eng. 2, 93-100 (1988)

        Tomii, K. and Kanehisa, M.;  Analysis of amino acid indices and
           mutation matrices for sequence comparison and structure
           prediction of proteins.  Protein Eng. 9, 27-36 (1996).

        Kawashima, S., Ogata, H., and Kanehisa, M.;  AAindex: amino acid
           index database.  Nucleic Acids Res. 27, 368-369 (1999).

        Kawashima, S. and Kanehisa, M.;  AAindex: amino acid index
           database.  Nucleic Acids Res. 28, 374 (2000).

        Last update: 2017/02/01
        """,
        [
            "https://www.genome.jp/ftp/db/community/aaindex/aaindex1",
            "https://www.genome.jp/ftp/db/community/aaindex/aaindex2",
            "https://www.genome.jp/ftp/db/community/aaindex/aaindex3",
        ],
        [
            "560f102dbee8ca9837186c4b2b93e253ea6efc58c88617234cef4e54bd61f676"
            "d3eb30c075452aa737b111a3d127358d0b669d0e847fca967b253a05c7796aa9"
            "8c9b9b0e7b7937b4674585b25f3f9370499f219b2583b1204fc9a0a5f06b1671"
        ]
    ))

    if !isfile(joinpath(datadep"AAindex", "aaindex.jld2"))
        preprocess_aaindex()
    end
end


text:
str = readline()
for c in str
    if islowercase(c)
        print(uppercase(c))
    else
        print(lowercase(c))
    end
end
println()


text:
using ReachabilityAnalysis
Dict([
    "opts" => Dict(:T=>20.0, :alg=>TMJets(abs_tol=1e-18, orderT=10, orderQ=4, max_steps=10_000)),
    "X0"   => Hyperrectangle(low=[0.19, 0.89], high=[0.2, 0.9]),
    "info" => Dict(
        "name" => "doublegyre",
        "dim" => 2,
        "linear" => false,
        "hybrid" => false
    )
])


text:
function _precompile_()
    ccall(:jl_generating_output, Cint, ()) == 1 || return nothing
    Base.precompile(Tuple{Core.kwftype(typeof(MatLang.freqspaceM)),NamedTuple{(:dim,),Tuple{Int64}},typeof(freqspaceM),Int64})
    Base.precompile(Tuple{Core.kwftype(typeof(MatLang.ndgridM)),NamedTuple{(:dim,),Tuple{Int64}},typeof(ndgridM),StepRange{Int64,Int64}})
    Base.precompile(Tuple{Core.kwftype(typeof(MatLang.randM)),NamedTuple{(:like,),Tuple{Array{Int64,2}}},typeof(randM),Int64})
    Base.precompile(Tuple{Core.kwftype(typeof(MatLang.sortperm)),NamedTuple{(:dims,),Tuple{Int64}},typeof(sortperm),Array{Int64,2}})
    Base.precompile(Tuple{typeof(blkdiagM),Array{Float64,2},Array{Float64,2},Vararg{Array{Float64,2},N} where N})
    Base.precompile(Tuple{typeof(catM),Int64,Array{Float64,2},Vararg{Array{Float64,2},N} where N})
    Base.precompile(Tuple{typeof(catM),Int64,Array{Float64,3},Vararg{Array{Float64,3},N} where N})
    Base.precompile(Tuple{typeof(diagM),Array{Int64,1}})
    Base.precompile(Tuple{typeof(eyeM),Int64})
    Base.precompile(Tuple{typeof(flipM),Array{Int64,1}})
    Base.precompile(Tuple{typeof(flipM),Array{Int64,2}})
    Base.precompile(Tuple{typeof(logspaceM),Int64,Int64,Int64})
    Base.precompile(Tuple{typeof(logspaceM),Int64,Int64,Symbol})
    Base.precompile(Tuple{typeof(logspaceM),Int64,Int64})
    Base.precompile(Tuple{typeof(meshgridM),StepRange{Int64,Int64}})
    Base.precompile(Tuple{typeof(meshgridM),UnitRange{Int64},UnitRange{Int64},UnitRange{Int64}})
    Base.precompile(Tuple{typeof(ndgridM),UnitRange{Int64},UnitRange{Int64},UnitRange{Int64}})
    Base.precompile(Tuple{typeof(randM),Symbol,Int64})
    Base.precompile(Tuple{typeof(randM),Type,Int64,Int64})
    Base.precompile(Tuple{typeof(repelemM),Array{Int64,2},Int64})
    Base.precompile(Tuple{typeof(repmatM),Array{Int64,1},Array{Int64,1}})
    Base.precompile(Tuple{typeof(repmatM),Int64,Int64,Int64})
    Base.precompile(Tuple{typeof(sizeM),Array{Float64,3},Int64,Int64})
    Base.precompile(Tuple{typeof(sizeM),Symbol,Array{Float64,3},Int64,Int64})
    Base.precompile(Tuple{typeof(sizeM),Symbol,Array{Int64,2}})
    Base.precompile(Tuple{typeof(sortM),Array{Int64,1}})
    Base.precompile(Tuple{typeof(sortM),Array{Int64,2}})
    Base.precompile(Tuple{typeof(sortM),Array{Integer,3},Int64})
    Base.precompile(Tuple{typeof(squeezeM),Array{Float64,3}})
    Base.precompile(Tuple{typeof(trueM),Array{Int64,1}})
    Base.precompile(Tuple{typeof(zerosM),Array{Int64,1}})
    Base.precompile(Tuple{typeof(zerosM),Type{Int32},Array{Int64,1}})
    isdefined(MatLang, Symbol("##sortperm#31")) && Base.precompile(Tuple{getfield(MatLang, Symbol("##sortperm#31")),Int64,Bool,typeof(sortperm),Array{Int64,2}})
end

max_stars_repo_name: juliamatlab/MatLang

max_issues_repo_name: juliamatlab/MatLang

max_forks_repo_name: juliamatlab/MatLang


text:
"""
Shorthand for `isnothing(optional) ? fallback : optional`
"""
@inline ifnothing(optional, fallback) = isnothing(optional) ? fallback : optional

function alignedbboxnode!(
    suggestedbbox::Node{BBox},
    computedsize::Node{NTuple{2, Optional{Float32}}},
    alignment::Node,
    sizeattrs::Node,
    autosizenode::Node{NTuple{2, Optional{Float32}}})

    finalbbox = Node(BBox(0, 100, 0, 100))

    onany(suggestedbbox, alignment, computedsize) do sbbox, al, csize

        bw = width(sbbox)
        bh = height(sbbox)

        # we only passively retrieve sizeattrs here because if they change
        # they also trigger computedsize, which triggers this node, too
        # we only need to know here if there are relative sizes given, because
        # those can only be computed knowing the suggestedbbox
        widthattr, heightattr = sizeattrs[]

        cwidth, cheight = csize

        w = if isnothing(cwidth)
            @match widthattr begin
                wa::Relative => wa.x * bw
                wa::Nothing => bw
                wa::Auto => if isnothing(autosizenode[][1])
                        # we have no autowidth available anyway
                        # take suggested width
                        bw
                    else
                        # use the width that was auto-computed
                        autosizenode[][1]
                    end
                wa => error("At this point, if computed width is not known,
                widthattr should be a Relative or Nothing, not $wa.")
            end
        else
            cwidth
        end

        h = if isnothing(cheight)
            @match heightattr begin
                ha::Relative => ha.x * bh
                ha::Nothing => bh
                ha::Auto => if isnothing(autosizenode[][2])
                        # we have no autoheight available anyway
                        # take suggested height
                        bh
                    else
                        # use the height that was auto-computed
                        autosizenode[][2]
                    end
                ha => error("At this point, if computed height is not known,
                heightattr should be a Relative or Nothing, not $ha.")
            end
        else
            cheight
        end

        # how much space is left in the bounding box
        rw = bw - w
        rh = bh - h

        xshift = @match al[1] begin
            :left => 0.0f0
            :center => 0.5f0 * rw
            :right => rw
            x::Real => x * rw
            x => error("Invalid horizontal alignment $x (only Real or :left, :center, or :right allowed).")
        end

        yshift = @match al[2] begin
            :bottom => 0.0f0
            :center => 0.5f0 * rh
            :top => rh
            x::Real => x * rh
            x => error("Invalid vertical alignment $x (only Real or :bottom, :center, or :top allowed).")
        end

        # align the final bounding box in the layout bounding box
        l = left(sbbox) + xshift
        b = bottom(sbbox) + yshift
        r = l + w
        t = b + h

        newbbox = BBox(l, r, b, t)
        # if finalbbox[] != newbbox
        #     finalbbox[] = newbbox
        # end
        finalbbox[] = newbbox
    end

    finalbbox
end

function computedsizenode!(sizeattrs, autosizenode::Node{NTuple{2, Optional{Float32}}})

    # set up csizenode with correct type manually
    csizenode = Node{NTuple{2, Optional{Float32}}}((nothing, nothing))

    onany(sizeattrs, autosizenode) do sizeattrs, autosize

        wattr, hattr = sizeattrs
        wauto, hauto = autosize

        wsize = computed_size(wattr, wauto)
        hsize = computed_size(hattr, hauto)

        csizenode[] = (wsize, hsize)
    end

    # trigger first value
    sizeattrs[] = sizeattrs[]

    csizenode
end

function computed_size(sizeattr, autosize)
    ms = @match sizeattr begin
        sa::Nothing => nothing
        sa::Real => sa
        sa::Fixed => sa.x
        sa::Relative => nothing
        sa::Auto => if sa.trydetermine
                # if trydetermine we report the autosize to the layout
                autosize
            else
                # but not if it's false, this allows for single span content
                # not to shrink its column or row, like a small legend next to an
                # axis or a super title over a single axis
                nothing
            end
        sa => error("""
            Invalid size attribute $sizeattr.
            Can only be Nothing, Fixed, Relative, Auto or Real""")
    end
end

function sizenode!(widthattr::Node, heightattr::Node)
    sizeattrs = Node{Tuple{Any, Any}}((widthattr[], heightattr[]))
    onany(widthattr, heightattr) do w, h
        sizeattrs[] = (w, h)
    end
    sizeattrs
end

function sceneareanode!(finalbbox, limits, aspect)

    scenearea = Node(IRect(0, 0, 100, 100))

    onany(finalbbox, limits, aspect) do bbox, limits, aspect

        w = width(bbox)
        h = height(bbox)
        # mw = min(w, maxsize[1])
        # mh = min(h, maxsize[2])
        # as = mw / mh
        as = w / h
        mw, mh = w, h

        if aspect isa AxisAspect
            aspect = aspect.aspect
        elseif aspect isa DataAspect
            aspect = limits.widths[1] / limits.widths[2]
        end

        if !isnothing(aspect)
            if as >= aspect
                # too wide
                mw *= aspect / as
            else
                # too high
                mh *= as / aspect
            end
        end

        restw = w - mw
        resth = h - mh

        # l = left(bbox) + alignment[1] * restw
        # b = bottom(bbox) + alignment[2] * resth
        l = left(bbox) + 0.5f0 * restw
        b = bottom(bbox) + 0.5f0 * resth

        newbbox = BBox(l, l + mw, b, b + mh)

        # only update scene if pixel positions change
        new_scenearea = IRect2D(newbbox)
        if new_scenearea != scenearea[]
            scenearea[] = new_scenearea
        end
    end

    scenearea
end

function create_suggested_bboxnode(n::Nothing)
    Node(BBox(0, 100, 0, 100))
end

function create_suggested_bboxnode(tup::Tuple)
    Node(BBox(tup...))
end

function create_suggested_bboxnode(bbox::AbstractPlotting.Rect2D)
    Node(BBox(bbox))
end

function create_suggested_bboxnode(node::Node{BBox})
    node
end

function roundedrectvertices(rect, cornerradius, cornersegments)
    cr = cornerradius
    csegs = cornersegments

    cr = min(width(rect) / 2, height(rect) / 2, cr)

    # inner corners
    ictl = topleft(rect) .+ Point2(cr, -cr)
    ictr = topright(rect) .+ Point2(-cr, -cr)
    icbl = bottomleft(rect) .+ Point2(cr, cr)
    icbr = bottomright(rect) .+ Point2(-cr, cr)

    # check if corners touch so we can remove one vertex that is doubled
    wtouching = width(rect) / 2 == cr
    htouching = height(rect) / 2 == cr

    cstr = if wtouching
        anglepoint.(Ref(ictr), LinRange(0, pi/2, csegs), cr)
    else
        anglepoint.(Ref(ictr), LinRange(0, pi/2, csegs)[1:end-1], cr)
    end
    cstl = if htouching
        anglepoint.(Ref(ictl), LinRange(pi/2, pi, csegs), cr)
    else
        anglepoint.(Ref(ictl), LinRange(pi/2, pi, csegs)[1:end-1], cr)
    end
    csbl = if wtouching
        anglepoint.(Ref(icbl), LinRange(pi, 3pi/2, csegs), cr)
    else
        anglepoint.(Ref(icbl), LinRange(pi, 3pi/2, csegs)[1:end-1], cr)
    end
    csbr = if htouching
        anglepoint.(Ref(icbr), LinRange(3pi/2, 2pi, csegs), cr)
    else
        anglepoint.(Ref(icbr), LinRange(3pi/2, 2pi, csegs)[1:end-1], cr)
    end
    arr = [cstr; cstl; csbl; csbr]
end

function anglepoint(center::Point2, angle::Real, radius::Real)
    Ref(center) .+ Ref(Point2(cos(angle), sin(angle))) .* radius
end

function enlarge(bbox::BBox, l, r, b, t)
    BBox(left(bbox) - l, right(bbox) + r, bottom(bbox) - b, top(bbox) + t)
end

function center(bbox::BBox)
    Point2f0((right(bbox) + left(bbox)) / 2, (top(bbox) + bottom(bbox)) / 2)
end

"""
Converts a point in fractions of rect dimensions into real coordinates.
"""
function fractionpoint(bbox::BBox, point::T) where T <: Point2
    T(left(bbox) + point[1] * width(bbox), bottom(bbox) + point[2] * height(bbox))
end

function tightlimits!(la::LAxis)
    la.xautolimitmargin = (0, 0)
    la.yautolimitmargin = (0, 0)
    autolimits!(la)
end

function tightlimits!(la::LAxis, sides::Union{Left, Right, Bottom, Top}...)
    for s in sides
        tightlimits!(la, s)
    end
end

function tightlimits!(la::LAxis, ::Left)
    la.xautolimitmargin = Base.setindex(la.xautolimitmargin[], 0.0, 1)
    autolimits!(la)
end

function tightlimits!(la::LAxis, ::Right)
    la.xautolimitmargin = Base.setindex(la.xautolimitmargin[], 0.0, 2)
    autolimits!(la)
end

function tightlimits!(la::LAxis, ::Bottom)
    la.yautolimitmargin = Base.setindex(la.yautolimitmargin[], 0.0, 1)
    autolimits!(la)
end

function tightlimits!(la::LAxis, ::Top)
    la.yautolimitmargin = Base.setindex(la.yautolimitmargin[], 0.0, 2)
    autolimits!(la)
end

function Base.foreach(f::Function, contenttype::Type, layout::GridLayout; recursive = true)
    for c in layout.content
        if recursive && c.content isa GridLayout
            foreach(f, contenttype, c.content)
        elseif c.content isa contenttype
            f(c.content)
        end
    end
end

"""
Swaps or rotates the layout positions of the given elements to their neighbor's.
"""
function swap!(layout_elements...)
    gridcontents = [le.layoutnodes.gridcontent for le in layout_elements]

    # copy relevant fields before gridcontents are mutated
    parents = map(gc -> gc.parent, gridcontents)
    spans = map(gc -> gc.span, gridcontents)
    sides = map(gc -> gc.side, gridcontents)

    for (gc, parent, span, side) in zip(circshift(gridcontents, 1), parents, spans, sides)
        parent[span.rows, span.cols, side] = gc.content
    end
end


text:
using StableRNGs
using MLJScikitLearnInterface
using Test
import MLJBase
import ScikitLearn
import MLJBase: target_scitype, input_scitype, output_scitype

const MB = MLJBase

include("testutils.jl")
include("macros.jl")

println("\nlinear-regressors");       include("models/linear-regressors.jl")
println("\nlinear-regressors-multi"); include("models/linear-regressors-multi.jl")

println("\nlinear-classifiers");      include("models/linear-classifiers.jl")
println("\ndiscriminant-analysis");   include("models/discriminant-analysis.jl")
println("\ngaussian-process");        include("models/gaussian-process.jl")
println("\nensemble");                include("models/ensemble.jl")
println("\nclustering");              include("models/clustering.jl")


text:
function calcarr(xa::Float64, p::Int64, start = 0.0, fin = 4.0)
  λarr = collect(linspace(start, fin, p))
  const iterconst = 50000
  const ϵ = 1e-6
  xarr = Array{Array{Float64,1},1}(p)
  yarr = Array{Array{Float64,1},1}(p)
  for j = 1:p
    λ = λarr[j]
    x = xa
    xarr[j] = Array{Float64,1}(0)
    yarr[j] = Array{Float64,1}(0)
    for i = 1:iterconst
      x = λ*x*(1-x)
    end
    push!(yarr[j], x)
    found = false
    while !found # while we still have new values
      x = λ*x*(1-x) # iterate once more
      for y in yarr[j] # iterate through all found values
        if almostequal(x, y, ϵ) # and check if we already have the new one
          found = true
          break
        end
      end
      if !found
        push!(yarr[j], x) # if not, add it to the array
      end
    end
    push!(xarr[j], λ)
  end
  return xarr, yarr
end

# function calcarrpmap(xa::Float64, p::Int64, start = 0.0, fin = 4.0)
#
# end
#
# xa = 0.5
# p = 1000
# start = 0.0
# fin = 4.0
#
# const chunks = 10
# step = div(p, 100)
# λarr = collect(linspace(start, fin, p))
# const iterconst = 50000
# const ϵ = 1e-6
# xarr = Array{Array{Float64, 1}, 1}(p)
# yarr = Array{Array{Float64, 1}, 1}(p)


text:
"""
    TMJets20{N, DM<:AbstractDisjointnessMethod} <: AbstractContinuousPost

Validated integration using Taylor models.

### Fields

- `orderQ`       -- (optional, default: `2`) order of the Taylor models for jet transport variables
- `orderT`       -- (optional, default: `8`) order of the Taylor model in time
- `abstol`       -- (optional, default: `1e-10`) absolute tolerance
- `maxsteps`     -- (optional, default: `2000`) maximum number of steps in the
                     validated integration ``x' = f(x)``
- `adaptive`     -- (optional, default: `true`) if `true`, try decreasing the absolute
                     tolerance each time step validation fails, until `min_abs_tol` is reached
- `minabstol`    -- (optional, default: `1e-29`) minimum absolute tolerance for the adaptive algorithm
- `disjointness` -- (optional, default: `ZonotopeEnclosure()`) defines the method to
                     perform the disjointness check between the taylor model flowpipe and the invariant

### Notes

The argument `disjointness` allows to control how are disjointness checks
computed, in the case where the invariant is not universal. In particular,
`ZonotopeEnclosure()` pre-processes the taylor model with a zonotopic
overapproximation, then performs the disjointness check with that zonotope and the
invariant. For other options, see the documentation of `AbstractDisjointnessMethod`.

This algorithm is an adaptation of the implementation in `TaylorModels.jl`
(see copyright license in the file `reach.jl` of the current folder). The package
`TaylorIntegration.jl` is used for jet-transport of ODEs using the Taylor method,
and `TaylorSeries.jl` is used to work with truncated Taylor series.
"""
@with_kw struct TMJets20{N, DM<:AbstractDisjointnessMethod} <: AbstractContinuousPost
    orderQ::Int=DEFAULT_ORDER_Q_TMJETS
    orderT::Int=DEFAULT_ORDER_T_TMJETS
    abstol::N=DEFAULT_ABS_TOL_TMJETS
    maxsteps::Int=DEFAULT_MAX_STEPS_TMJETS
    adaptive::Bool=true
    minabstol::N=1e-29
    disjointness::DM=ZonotopeEnclosure()
end

const DEFAULT_MAX_STEPS_TMJETS = 2000
const DEFAULT_ABS_TOL_TMJETS = 1e-10
const DEFAULT_ORDER_T_TMJETS = 8
const DEFAULT_ORDER_Q_TMJETS = 2

numtype(::TMJets20{N}) where {N} = N
rsetrep(::TMJets20{N}) where {N} = TaylorModelReachSet{N}

include("post.jl")
include("reach.jl")


text:
using PictureSegmentation
using Test

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


text:
# additional tests for config

foofig(p, s) = (write(joinpath(p, "src", "config.md"), s); J.process_config())

@testset "config" begin
    p = joinpath(D, "..", "__tst_config")
    isdir(p) && rm(p; recursive=true, force=true)
    mkdir(p); cd(p);
    J.FOLDER_PATH[] = pwd()
    J.set_paths!()
    mkdir(joinpath(p, "src"))
    # ================================
    # asssignments go to GLOBAL
    foofig(p, raw"""
    @def var = 5
    """)
    @test haskey(J.GLOBAL_PAGE_VARS, "var")
    @test J.GLOBAL_PAGE_VARS["var"][1] == 5

    # lxdefs go to GLOBAL
    foofig(p, raw"""
    \newcommand{\hello}{goodbye}
    """)
    @test haskey(J.GLOBAL_LXDEFS, "\\hello")
    @test J.GLOBAL_LXDEFS["\\hello"].def == "goodbye"

    # combination of lxdefs
    foofig(p, raw"""
    \newcommand{\hello}{goodbye}
    \newcommand{\hellob}{\hello}
    \newcommand{\helloc}{\hellob}
    """)

    @test J.GLOBAL_LXDEFS["\\hello"].from < J.GLOBAL_LXDEFS["\\hello"].to <
          J.GLOBAL_LXDEFS["\\hellob"].from < J.GLOBAL_LXDEFS["\\hellob"].to <
          J.GLOBAL_LXDEFS["\\helloc"].from < J.GLOBAL_LXDEFS["\\helloc"].to

    @test jd2html(raw"""\helloc"""; dir=p, internal=true) == "goodbye"

    # ================================
    # go back and cleanup
    cd(R); rm(p; recursive=true, force=true)
end


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

JLLWrappers.@generate_main_file_header("HELICS")
JLLWrappers.@generate_main_file("HELICS", UUID("ef3b0bb0-9dc6-5204-90f3-946fd7d0da3e"))
end  # module HELICS_jll


text:
"""
Add a string to the beginning of each element of each key of a dict
"""
function prepend!(d::Dict,pre::String)
  dnew = deepcopy(d)
  for key in keys(d)
    for (i,dname) in enumerate(d[key])
      dnew[key][i] = pre*dname
    end
  end
  return dnew
end

function subs(x::Array{Sym},subin::Array{Sym},subout::Array{Sym})
  y = deepcopy(x)
  for i in eachindex(subin)
    y = subs(y,subin[i],subout[i])
  end
  return y
end

function subs(x::Sym,subin::Array{Sym},subout::Array{Sym})
  y = deepcopy(x)
  for i in eachindex(subin)
    y = subs(y,subin[i],subout[i])
  end
  return y
end

function *(x::Array{String},val::String)
  for i in eachindex(x)
    x[i] = x[i]*val
  end
  return x
end

function *(val::String,x::Array{String})
  for i in eachindex(x)
    x[i] = val*x[i]
  end
  return x
end

# function convert(::Type{Array{Sym,1}},x::Sym)
#   x = [Sym(x)]
# end


text:
# Autogenerated wrapper script for Prjtrellis_jll for armv7l-linux-gnueabihf-cxx11
export ecpbram, ecpmulti, ecpmulti, ecppack, ecppll, ecpunpack, libtrellis, pytrellis

using Python_jll
using boost_jll
JLLWrappers.@generate_wrapper_header("Prjtrellis")
JLLWrappers.@declare_executable_product(ecpbram)
JLLWrappers.@declare_executable_product(ecpmulti)
JLLWrappers.@declare_executable_product(ecpmulti)
JLLWrappers.@declare_executable_product(ecppack)
JLLWrappers.@declare_executable_product(ecppll)
JLLWrappers.@declare_executable_product(ecpunpack)
JLLWrappers.@declare_library_product(libtrellis, "libtrellis.so")
JLLWrappers.@declare_library_product(pytrellis, "pytrellis.so")
function __init__()
    JLLWrappers.@generate_init_header(Python_jll, boost_jll)
    JLLWrappers.@init_executable_product(
        ecpbram,
        "bin/ecpbram",
    )

    JLLWrappers.@init_executable_product(
        ecpmulti,
        "bin/ecpmulti",
    )

    JLLWrappers.@init_executable_product(
        ecpmulti,
        "bin/ecpmulti",
    )

    JLLWrappers.@init_executable_product(
        ecppack,
        "bin/ecppack",
    )

    JLLWrappers.@init_executable_product(
        ecppll,
        "bin/ecppll",
    )

    JLLWrappers.@init_executable_product(
        ecpunpack,
        "bin/ecpunpack",
    )

    JLLWrappers.@init_library_product(
        libtrellis,
        "lib/trellis/libtrellis.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        pytrellis,
        "lib/trellis/pytrellis.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
"""
    getExternalCycle(model::CAGD.Model, dim::Int, cells::Array{Int,1})::Int

Evaluates wich `dim`-coboundary between `cells` is the external one.

The function expects `cells` as the indices of ALL the `dim`-coboundary
within the same component and evaluate wich of those is the external one.

It returns the index of the external cell within the `cells`.
"""
function getExternalCycle(model::CAGD.Model, dim::Int, cells::Array{Int,1})::Int

    return 1
end

"""
    buildContainmentGraph(
        model::CAGD.Model,
        dim::Int,
        externals::Array{Int,1},
        internals::Array{Array{Int,1},1}
    )::SparseMatrixCSC{Int, Int}

Evaluates the strict containment relation between `dim`-cells.

The method expects `externals[i]` and `internals[i]` as indices respectivelly
of the external and internal cycles of the `i`-th component out of `ne`.

Returns a `ne × ne` sparse `graph` where `graph[i, j] = k` if `i`-th component
includes `j`-th one within the `k`-th cycle (internal cycle of `i`-th component)
"""
function buildContainmentGraph(model, dim, externals, internals)
    ne = length(externals)

    # Bounding box for all the `dim`-cells
    bboxes = CAGD.getModelBoundingBoxes(model, dim)

    # graph[i, j] = 1 if bounding box of external[i] contains the bounding box of external[j]
    graph = SparseArrays.sparse([], [], Bool.([]), ne, ne)
    for c = 1 : ne  for r = 1 : ne
        if CAGD.isBoundingBoxIncluded(bboxes[externals[c]], bboxes[externals[r]])
            graph[r, c] = true
        end
    end  end

    # define a point per contained cell to be tested out
    pts = zeros(3, ne)
    containers, contained, _ = SparseArrays.findnz(graph)
    for c in unique(contained)
        pts[:, c] = CAGD.getModelCellVertices(model, dim, externals[c])[1]
    end

    # graph[r, c] = 1 if a point of c-th (external) is in r-th cell (external)
    for (r, c) in zip(containers, contained)
        if CAGD.pointInCell(model, dim, pts[c], externals[r], signed = false) >= 0
            SparseArrays.dropstored!(graph, r, c)
        end
    end

    # Transitive Reduction of graph
    containers, contained, _ = SparseArrays.findnz(graph)
    for (out, mid) in zip(containers, contained)
        for inn in SparseArrays.nzind(graph[mid, :])
            SparseArrays.dropstored!(graph, out, inn);
        end
    end

    # graph[r, c] = k if internal[r][k] contains external[c]
    containers, contained, _ = SparseArrays.findnz(graph)
    for (container, contained) in zip(containers, contained)
        for doContain in internals[container]
            if CAGD.isBoundingBoxIncluded(bboxes[externals[contained]], bboxes[doContain]) &
                    (CAGD.pointInCell(model, dim, pts[contained], doContain, signed = false) >= 0)
                graph[container, contained] = doContain
                break
            end
        end
    end

    return graph
end

"""
    assembleModel!(
        model::CAGD.Model,
        dim::Int,
        externals::Array{Int,1},
        containment_graph::SparseMatrixCSC{Int, Int},
        doExternal = false
    )::Nothing

Performs the adjoining of the outer cycles of inner components to the outer ones

Assemble on `model` the `ne` `dim`-cycles `externals` on inner cycle specified
by `contained_graph`. The latter is expected to be a `ne × ne` `SparseMatrixCSC`
where `containment_graph[i, j] = k` if `i`-th component contains the `j`-th
component within cycle `k`.

If `doExternal` is set to true, then the outer cycle is also adjoint.
"""
function assembleModel!(model, dim, externals, containment_graph, doExternal = false)
    containers, contained, containing = SparseArrays.findnz(containment_graph)

    # Adjoining outer cell of inner components to outer components
    for (out, inn) in zip(containing, contained)
        model.T[dim][out, :] += model.T[dim][externals[inn], :]
    end

    # Retrieve the external cycle (if needed)
    if doExternal
        external = SparseArrays.spzeros(Int8, size(model, dim, 2))
        for out in setdiff(1 : length(externals), unique(containers))
            external += model.T[dim][externals[out], :]
        end
        CAGD.addModelCells!(model, dim, external)
    end

    # Removing outer cells of components
    CAGD.deleteModelCells!(model, dim, externals)
    return
end


text:
function Base.:^(a::UncertainValueDataset, b::UncertainValueDataset; n = 30000)
    N = length(a)
    n_vals_b = length(b)

    if N != n_vals_b
        throw(ArgumentError("Dataset lengths do not match ($N, $n_vals_b)"))
    end

    UncertainValueDataset([a[i] ^ b[i] for i = 1:N])
end

function Base.:^(a::T, b::UncertainValueDataset; n = 30000) where {T <: Number}
    N = length(a)
    UncertainValueDataset([a ^ b[i] for i = 1:N])
end

function Base.:^(a::Vector{T}, b::UncertainValueDataset; n = 30000) where {T <: Number}
    N = length(a)
    n_vals_b = length(b)

    if N != n_vals_b
        throw(ArgumentError("Dataset lengths do not match ($N, $n_vals_b)"))
    end

    UncertainValueDataset([a[i] ^ b[i] for i = 1:N])
end

function Base.:^(a::UncertainValueDataset, b::T; n = 30000) where {T <: Number}
    N = length(a)

    UncertainValueDataset([a[i] ^ b for i = 1:N])
end

function Base.:^(a::UncertainValueDataset, b::Vector{T}; n = 30000) where {T <: Number}
    N = length(a)
    n_vals_b = length(b)

    if N != n_vals_b
        throw(ArgumentError("Dataset lengths do not match ($N, $n_vals_b)"))
    end

    UncertainValueDataset([a[i] ^ b[i] for i = 1:N])
end


text:
using Viznet
using Viznet: inner_most_containers, put_edge!, put_node!, empty_cache!, nedge, nnode,
    EDGE_CACHE, NODE_CACHE
using Test
using Compose

@testset "basic" begin
    empty_cache!()
    c = compose(context(), compose(context(), line()))
    n = compose(context(), compose(context(), circle()))
    ics = []
    inner_most_containers(x->!isempty(x.form_children) && push!(ics, x), c)
    @test length(ics) == 1
    ic = ics[]
    @test first(ic.form_children) isa Compose.Form{<:Compose.LinePrimitive}
    put_edge!(c, (0.3, 0.4), (2.3, 1.4))
    put_edge!(c, (0.3, 0.5), (2.3, 1.4))
    put_node!(n, (0.3, 0.5))
    @test nnode() == 1
    @test nedge() == 2
    empty_cache!()
    @test nnode() == 0
    @test nedge() == 0
end

@testset "flush edges" begin
    empty_cache!()
    c = compose(context(), compose(context(), line()))
    c >> ((0.3, 0.4), (2.3, 1.4))
    c >> ((0.3, 0.8), (2.3, 1.9))
    c0 = bondstyle(:default)
    c0 >> ((0.3, 0.8), (2.3, 2.9))
    lst = flush!(EDGE_CACHE)
    @test length(lst) == 2
end

@testset "flush nodes" begin
    empty_cache!()
    c = compose(context(), compose(context(), circle(0.2, 0.2, 2.0)))
    c >> (0.3, 0.4)
    c >> (0.3, 0.8)
    c0 = nodestyle(:default)
    c0 >> (0.3, 0.8)
    lst = flush!(NODE_CACHE)
    @test length(lst) == 2
end

@testset "system test" begin
    empty_cache!()
    nb = compose(context(), stroke("red"), circle(0.0, 0.0, 0.02))
    eb = compose(context(), stroke("green"), line())

    x = nb >> (0.2, 0.3)
    y = nb >> (0.6, 0.6)
    eb >> (x, y)
    @test flush!() isa Context
end

@testset "similar shape" begin
    for p2 in [
        polygon([(0.1, 0.2), (0.2, 0.3), (0.3,0.3)]),
        circle(0.1, 0.2, 0.3),
        rectangle(0.1,0.1, 0.2, 0.2),
        line(),
        arc(0.1, 0.1, 0.2, π, π/2),
        curve((0.0, 0.1), (0.2, 0.2), (0.3, 0.4), (0.5, 0.0)),
        ]
        nodes = Viznet.similar(p2, [(0.1,0.1), (0.2,0.2), (0.3,0.3), (0.4, 0.4)])
        @test nodes.primitives |> length == 4
    end
    for l2 in [
        line(),
        arc(0.1, 0.1, 0.2, π, π/2),
        curve((0.0, 0.1), (0.2, 0.2), (0.3, 0.4), (0.5, 0.0)),
        ]
        lines = Viznet.similar(l2, [((0.1,0.1), (0.2,0.2)), ((0.3,0.3), (0.4, 0.4))])
        @test lines.primitives |> length == 2
    end
end


text:
using Test, Pluto

function testnotebook(input)
    # de-markdown eventual cells with Pkg.develop and write
    # to pluto-tmp.jl
    notebook=Pluto.load_notebook_nobackup(input)
    pkgcellfortest=findfirst(c->occursin("Pkg.develop",c.code),notebook.cells)
    if  pkgcellfortest!=nothing
        # de-markdown pkg cell
        notebook.cells[pkgcellfortest].code=replace(notebook.cells[pkgcellfortest].code,"md"=>"")
        notebook.cells[pkgcellfortest].code=replace(notebook.cells[pkgcellfortest].code,"\"\"\""=>"")
        notebook.cells[pkgcellfortest].code=replace(notebook.cells[pkgcellfortest].code,";"=>"")
        @info "Pkg cell: $(pkgcellfortest)\n$(notebook.cells[pkgcellfortest].code)"
        Pluto.save_notebook(notebook,"pluto-tmp.jl")
        input="pluto-tmp.jl"
        sleep(1)
    end

    # run notebook and check for cell errors
    session = Pluto.ServerSession();
    notebook = Pluto.SessionActions.open(session, input; run_async=false)
    errored=false
    for c in notebook.cells
        if c.errored
            errored=true
            @error "Error in  $(c.cell_id): $(c.output.body[:msg])\n $(c.code)"
        end
    end
    !errored
end

@testset "notebooks" begin
    notebooks=["vtktest.jl","plotlytest.jl","plutovista.jl"]
    for notebook in notebooks
        input=joinpath(@__DIR__,"..","examples",notebook)
        @info "notebook: $(input)"
        @test testnotebook(input)
    end
end


text:
# syntax: proto3
using ProtoBuf
import ProtoBuf.meta

mutable struct ResourceHandleProto_DtypeAndShape <: ProtoType
    __protobuf_jl_internal_meta::ProtoMeta
    __protobuf_jl_internal_values::Dict{Symbol,Any}
    __protobuf_jl_internal_defaultset::Set{Symbol}

    function ResourceHandleProto_DtypeAndShape(; kwargs...)
        obj = new(meta(ResourceHandleProto_DtypeAndShape), Dict{Symbol,Any}(), Set{Symbol}())
        values = obj.__protobuf_jl_internal_values
        symdict = obj.__protobuf_jl_internal_meta.symdict
        for nv in kwargs
            fldname, fldval = nv
            fldtype = symdict[fldname].jtyp
            (fldname in keys(symdict)) || error(string(typeof(obj), " has no field with name ", fldname))
            values[fldname] = isa(fldval, fldtype) ? fldval : convert(fldtype, fldval)
        end
        obj
    end
end # mutable struct ResourceHandleProto_DtypeAndShape
const __meta_ResourceHandleProto_DtypeAndShape = Ref{ProtoMeta}()
function meta(::Type{ResourceHandleProto_DtypeAndShape})
    ProtoBuf.metalock() do
        if !isassigned(__meta_ResourceHandleProto_DtypeAndShape)
            __meta_ResourceHandleProto_DtypeAndShape[] = target = ProtoMeta(ResourceHandleProto_DtypeAndShape)
            allflds = Pair{Symbol,Union{Type,String}}[:dtype => Int32, :shape => TensorShapeProto]
            meta(target, ResourceHandleProto_DtypeAndShape, allflds, ProtoBuf.DEF_REQ, ProtoBuf.DEF_FNUM, ProtoBuf.DEF_VAL, ProtoBuf.DEF_PACK, ProtoBuf.DEF_WTYPES, ProtoBuf.DEF_ONEOFS, ProtoBuf.DEF_ONEOF_NAMES)
        end
        __meta_ResourceHandleProto_DtypeAndShape[]
    end
end
function Base.getproperty(obj::ResourceHandleProto_DtypeAndShape, name::Symbol)
    if name === :dtype
        return (obj.__protobuf_jl_internal_values[name])::Int32
    elseif name === :shape
        return (obj.__protobuf_jl_internal_values[name])::TensorShapeProto
    else
        getfield(obj, name)
    end
end

mutable struct ResourceHandleProto <: ProtoType
    __protobuf_jl_internal_meta::ProtoMeta
    __protobuf_jl_internal_values::Dict{Symbol,Any}
    __protobuf_jl_internal_defaultset::Set{Symbol}

    function ResourceHandleProto(; kwargs...)
        obj = new(meta(ResourceHandleProto), Dict{Symbol,Any}(), Set{Symbol}())
        values = obj.__protobuf_jl_internal_values
        symdict = obj.__protobuf_jl_internal_meta.symdict
        for nv in kwargs
            fldname, fldval = nv
            fldtype = symdict[fldname].jtyp
            (fldname in keys(symdict)) || error(string(typeof(obj), " has no field with name ", fldname))
            values[fldname] = isa(fldval, fldtype) ? fldval : convert(fldtype, fldval)
        end
        obj
    end
end # mutable struct ResourceHandleProto
const __meta_ResourceHandleProto = Ref{ProtoMeta}()
function meta(::Type{ResourceHandleProto})
    ProtoBuf.metalock() do
        if !isassigned(__meta_ResourceHandleProto)
            __meta_ResourceHandleProto[] = target = ProtoMeta(ResourceHandleProto)
            allflds = Pair{Symbol,Union{Type,String}}[:device => AbstractString, :container => AbstractString, :name => AbstractString, :hash_code => UInt64, :maybe_type_name => AbstractString, :dtypes_and_shapes => Base.Vector{ResourceHandleProto_DtypeAndShape}]
            meta(target, ResourceHandleProto, allflds, ProtoBuf.DEF_REQ, ProtoBuf.DEF_FNUM, ProtoBuf.DEF_VAL, ProtoBuf.DEF_PACK, ProtoBuf.DEF_WTYPES, ProtoBuf.DEF_ONEOFS, ProtoBuf.DEF_ONEOF_NAMES)
        end
        __meta_ResourceHandleProto[]
    end
end
function Base.getproperty(obj::ResourceHandleProto, name::Symbol)
    if name === :device
        return (obj.__protobuf_jl_internal_values[name])::AbstractString
    elseif name === :container
        return (obj.__protobuf_jl_internal_values[name])::AbstractString
    elseif name === :name
        return (obj.__protobuf_jl_internal_values[name])::AbstractString
    elseif name === :hash_code
        return (obj.__protobuf_jl_internal_values[name])::UInt64
    elseif name === :maybe_type_name
        return (obj.__protobuf_jl_internal_values[name])::AbstractString
    elseif name === :dtypes_and_shapes
        return (obj.__protobuf_jl_internal_values[name])::Base.Vector{ResourceHandleProto_DtypeAndShape}
    else
        getfield(obj, name)
    end
end

export ResourceHandleProto_DtypeAndShape, ResourceHandleProto


text:
using Kinetic, Plots, LinearAlgebra, JLD2, Flux
using KitBase.ProgressMeter: @showprogress
using Flux: onecold
pyplot()

cd(@__DIR__)
include("tools.jl")

function ev!(
    KS::SolverSet,
    ctr::T1,
    a1face::T2,
    a2face::T2,
    dt;
    mode = KitBase.symbolize(KS.set.flux)::Symbol,
    bc = KitBase.symbolize(KS.set.boundary),
) where {T1<:AbstractArray{ControlVolume2D2F,2},T2<:AbstractArray{Interface2D2F,2}}

    nx, ny, dx, dy = begin
        if KS.ps isa CSpace2D
            KS.ps.nr, KS.ps.nθ, KS.ps.dr, KS.ps.darc
        else
            KS.ps.nx, KS.ps.ny, KS.ps.dx, KS.ps.dy
        end
    end

    if firstindex(KS.pSpace.x[:, 1]) < 1
        idx0 = 1
        idx1 = nx + 1
    else
        idx0 = 2
        idx1 = nx
    end
    if firstindex(KS.pSpace.y[1, :]) < 1
        idy0 = 1
        idy1 = ny + 1
    else
        idy0 = 2
        idy1 = ny
    end

    # x direction
    @inbounds Threads.@threads for j = 1:ny
        for i = idx0:idx1
            vn = KS.vSpace.u .* a1face[i, j].n[1] .+ KS.vSpace.v .* a1face[i, j].n[2]
            vt = KS.vSpace.v .* a1face[i, j].n[1] .- KS.vSpace.u .* a1face[i, j].n[2]

            w = (ctr[i-1, j].w + ctr[i, j].w) ./ 2
            f = (ctr[i-1, j].h + ctr[i, j].h) ./ 2
            swx = (ctr[i, j].w - ctr[i-1, j].w) / (ks.ps.x[i, j] - ks.ps.x[i-1, j])
            swy = (ctr[i, j].w - ctr[i-1, j].w) / (ks.ps.y[i, j] - ks.ps.y[i-1, j])
            #sw = sqrt.(sw1.^2 + sw2.^2)
            prim = conserve_prim(w, ks.gas.γ)
            #τ = vhs_collision_time(prim, ks.gas.μᵣ, ks.gas.ω)
            regime = judge_regime(ks, f, prim, swx, swy)

            if regime == 2
                flux_kfvs!(
                    a1face[i, j].fw,
                    a1face[i, j].fh,
                    a1face[i, j].fb,
                    ctr[i-1, j].h .+ 0.5 .* dx[i-1, j] .* ctr[i-1, j].sh[:, :, 1],
                    ctr[i-1, j].b .+ 0.5 .* dx[i-1, j] .* ctr[i-1, j].sb[:, :, 1],
                    ctr[i, j].h .- 0.5 .* dx[i, j] .* ctr[i, j].sh[:, :, 1],
                    ctr[i, j].b .- 0.5 .* dx[i, j] .* ctr[i, j].sb[:, :, 1],
                    vn,
                    vt,
                    KS.vSpace.weights,
                    dt,
                    a1face[i, j].len,
                    ctr[i-1, j].sh[:, :, 1],
                    ctr[i-1, j].sb[:, :, 1],
                    ctr[i, j].sh[:, :, 1],
                    ctr[i, j].sb[:, :, 1],
                )
            else
                wL = local_frame(ctr[i-1, j].w, a1face[i, j].n[1], a1face[i, j].n[2])
                wR = local_frame(ctr[i, j].w, a1face[i, j].n[1], a1face[i, j].n[2])
                swL = local_frame(ctr[i-1, j].sw[:, 1], a1face[i, j].n[1], a1face[i, j].n[2])
                swR = local_frame(ctr[i, j].sw[:, 1], a1face[i, j].n[1], a1face[i, j].n[2])

                flux_gks!(
                    a1face[i, j].fw,
                    a1face[i, j].fh,
                    a1face[i, j].fb,
                    wL .+ 0.5 .* dx[i, j] .* swL,
                    wR .- 0.5 .* dx[i, j] .* swR,
                    vn,
                    vt,
                    ks.gas.K,
                    ks.gas.γ,
                    ks.gas.μᵣ,
                    ks.gas.ω,
                    dt,
                    0.5 .* dx[i-1, j],
                    0.5 .* dx[i, j],
                    a1face[i, j].len,
                    swL,
                    swR,
                )
            end

            a1face[i, j].fw .=
                global_frame(a1face[i, j].fw, a1face[i, j].n[1], a1face[i, j].n[2])
        end
    end

    # y direction
    @inbounds Threads.@threads for j = idy0:idy1
        for i = 1:nx
            vn = KS.vSpace.u .* a2face[i, j].n[1] .+ KS.vSpace.v .* a2face[i, j].n[2]
            vt = KS.vSpace.v .* a2face[i, j].n[1] .- KS.vSpace.u .* a2face[i, j].n[2]

            w = (ctr[i, j-1].w + ctr[i, j].w) ./ 2
            f = (ctr[i, j-1].h + ctr[i, j].h) ./ 2
            swx = (ctr[i, j].w - ctr[i, j-1].w) / (ks.ps.x[i, j] - ks.ps.x[i, j-1])
            swy = (ctr[i, j].w - ctr[i, j-1].w) / (ks.ps.y[i, j] - ks.ps.y[i, j-1])
            #sw = sqrt.(sw1.^2 + sw2.^2)
            prim = conserve_prim(w, ks.gas.γ)
            #τ = vhs_collision_time(prim, ks.gas.μᵣ, ks.gas.ω)
            regime = judge_regime(ks, f, prim, swx, swy)

            if regime == 2
                flux_kfvs!(
                    a2face[i, j].fw,
                    a2face[i, j].fh,
                    a2face[i, j].fb,
                    ctr[i, j-1].h .+ 0.5 .* dy[i, j-1] .* ctr[i, j-1].sh[:, :, 2],
                    ctr[i, j-1].b .+ 0.5 .* dy[i, j-1] .* ctr[i, j-1].sb[:, :, 2],
                    ctr[i, j].h .- 0.5 .* dy[i, j] .* ctr[i, j].sh[:, :, 2],
                    ctr[i, j].b .- 0.5 .* dy[i, j] .* ctr[i, j].sb[:, :, 2],
                    vn,
                    vt,
                    KS.vSpace.weights,
                    dt,
                    a2face[i, j].len,
                    ctr[i, j-1].sh[:, :, 2],
                    ctr[i, j-1].sb[:, :, 2],
                    ctr[i, j].sh[:, :, 2],
                    ctr[i, j].sb[:, :, 2],
                )
            else
                wL = local_frame(ctr[i, j-1].w, a2face[i, j].n[1], a2face[i, j].n[2])
                wR = local_frame(ctr[i, j].w, a2face[i, j].n[1], a2face[i, j].n[2])
                swL = local_frame(ctr[i, j-1].sw[:, 2], a2face[i, j].n[1], a2face[i, j].n[2])
                swR = local_frame(ctr[i, j].sw[:, 2], a2face[i, j].n[1], a2face[i, j].n[2])

                flux_gks!(
                    a2face[i, j].fw,
                    a2face[i, j].fh,
                    a2face[i, j].fb,
                    wL .+ 0.5 .* dy[i, j-1] .* swL,
                    wR .- 0.5 .* dy[i, j] .* swR,
                    vn,
                    vt,
                    ks.gas.K,
                    ks.gas.γ,
                    ks.gas.μᵣ,
                    ks.gas.ω,
                    dt,
                    0.5 .* dy[i, j-1],
                    0.5 .* dy[i, j],
                    a2face[i, j].len,
                    swL,
                    swR,
                )
            end

            a2face[i, j].fw .=
                global_frame(a2face[i, j].fw, a2face[i, j].n[1], a2face[i, j].n[2])
        end
    end

    KitBase.evolve_boundary!(KS, ctr, a1face, a2face, dt, mode, bc)

    return nothing

end

begin
    set = Setup(
        case = "cylinder",
        space = "2d2f2v",
        boundary = ["maxwell", "extra", "mirror", "mirror"],
        limiter = "minmod",
        cfl = 0.5,
        maxTime = 10.0, # time
    )
    ps = CSpace2D(1.0, 6.0, 30, 0.0, π, 50, 1, 1)
    vs = VSpace2D(-10.0, 10.0, 48, -10.0, 10.0, 48)
    gas = Gas(Kn = 1e-2, Ma = 5.0, K = 1.0)

    prim0 = [1.0, 0.0, 0.0, 1.0]
    prim1 = [1.0, gas.Ma * sound_speed(1.0, gas.γ), 0.0, 1.0]
    fw = (args...) -> prim_conserve(prim1, gas.γ)
    ff = function(args...)
        prim = conserve_prim(fw(args...), gas.γ)
        h = maxwellian(vs.u, vs.v, prim)
        b = h .* gas.K / 2 / prim[end]
        return h, b
    end
    bc = function(x, y)
        if abs(x^2 + y^2 - 1) < 1e-3
            return prim0
        else
            return prim1
        end
    end
    ib = IB2F(fw, ff, bc)

    ks = SolverSet(set, ps, vs, gas, ib)
end

ctr, a1face, a2face = init_fvm(ks, ks.ps)

t = 0.0
dt = timestep(ks, ctr, 0.0)
nt = ks.set.maxTime ÷ dt |> Int
res = zeros(4)
@showprogress for iter = 1:100#nt
    #rc!(ks, ctr)
    ev!(ks, ctr, a1face, a2face, dt)
    #evolve!(ks, ctr, a1face, a2face, dt)
    update!(ks, ctr, a1face, a2face, dt, res)

    for j = ks.ps.nθ÷2+1:ks.ps.nθ
        ctr[ks.ps.nr+1, j].w .= ks.ib.fw(6, 0)
        ctr[ks.ps.nr+1, j].prim .= conserve_prim(ctr[ks.ps.nr+1, j].w, ks.gas.γ)
        ctr[ks.ps.nr+1, j].sw .= 0.0
        ctr[ks.ps.nr+1, j].h .= maxwellian(ks.vs.u, ks.vs.v, ctr[ks.ps.nr+1, j].prim)
        ctr[ks.ps.nr+1, j].b .= ctr[ks.ps.nr+1, j].h .* ks.gas.K ./ 2 ./ ctr[ks.ps.nr+1, j].prim[end]
    end

    global t += dt

    if maximum(res) < 1e-6
        break
    end
end

begin
    sol = zeros(ks.ps.nr, ks.ps.nθ, 4)
    for i in axes(sol, 1), j in axes(sol, 2)
        sol[i, j, :] .= ctr[i, j].prim
        sol[i, j, end] = 1 / sol[i, j, end]
    end
    contourf(
        ps.x[1:ks.ps.nr, 1:ks.ps.nθ],
        ps.y[1:ks.ps.nr, 1:ks.ps.nθ],
        sol[:, :, 4],
        ratio = 1,
    )
end

max_stars_repo_name: CSMMLab/Flowmachine

max_issues_repo_name: CSMMLab/Flowmachine

max_forks_repo_name: CSMMLab/Flowmachine


text:
using RungeKutta
using Documenter
using DocumenterCitations
using Weave

bib = CitationBibliography("RungeKutta.bib")

module RungeKuttaWeaves
    using Markdown
    using Markdown: MD, Paragraph, LineBreak
    using Plots
    using PrettyTables
    using RungeKutta
    using RungeKutta.Tableaus: get_gauss_nodes, get_gauss_weights, get_gauss_coefficients
    using RungeKutta.Tableaus: get_lobatto_nodes, get_lobatto_weights,
                               get_lobatto_a_coefficients, get_lobatto_b_coefficients,
                               get_lobatto_c_coefficients, get_lobatto_c̄_coefficients,
                               get_lobatto_d_coefficients, get_lobatto_e_coefficients,
                               get_lobatto_f_coefficients, get_lobatto_g_coefficients
    using RungeKutta.Tableaus: get_radau_1_nodes, get_radau_1_weights, get_radau_1_coefficients,
                               get_radau_2_nodes, get_radau_2_weights, get_radau_2_coefficients
    using SymPy
    using SymPy: latex

    "Markdown-print Runge-Kutta tableau with SymPy coefficients."
    function Base.show(io::IO, ::MIME"text/markdown", tab::Tableau{Sym})
        show(io, "text/markdown", Markdown.parse("Runge-Kutta Tableau $(tab.name) with $(tab.s) stages and order $(tab.o):"))

        tab_arr = simplify.(convert(Matrix, tab))
        str_arr = latex.(tab_arr)
        str_arr[tab.s+1,1] = ""

        strio = IOBuffer()
        pretty_table(strio, str_arr,
                        backend = :latex,
                        vlines = [1],
                        hlines = [tab.s],
                        noheader = true)
        tab_latex = String(take!(strio))

        tab_markdown = replace(tab_latex, "tabular" => "array")
        tab_markdown = replace(tab_markdown, "\\begin{table}" => "```math")
        tab_markdown = replace(tab_markdown, "\\end{table}" => "```")

        print(io, tab_markdown * "\n")
    end
end

weave("src/gauss.jmd",
         out_path = "src",
         doctype = "github",
         mod = RungeKuttaWeaves)

weave("src/radauIA.jmd",
         out_path = "src",
         doctype = "github",
         mod = RungeKuttaWeaves)

weave("src/radauIIA.jmd",
         out_path = "src",
         doctype = "github",
         mod = RungeKuttaWeaves)

weave("src/lobatto.jmd",
         out_path = "src",
         doctype = "github",
         mod = RungeKuttaWeaves)

makedocs(bib;
    modules=[RungeKutta],
    authors="Michael Kraus",
    repo="https://github.com/JuliaGNI/RungeKutta.jl/blob/{commit}{path}#L{line}",
    sitename="RungeKutta.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://juliagni.github.io/RungeKutta.jl",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
        "Runge-Kutta Methods" => "rungekutta.md",
        "Tableau" => "tableau.md",
        "Diagnostics" => "diagnostics.md",
        "Gauß Methods" => "gauss.md",
        "Radau IA Methods" => "radauIA.md",
        "Radau IIA Methods" => "radauIIA.md",
        "Lobatto III Methods" => "lobatto.md",
        "Library" => "library.md",
        "Bibliography" => "bibliography.md",
    ],
)

deploydocs(;
    repo="github.com/JuliaGNI/RungeKutta.jl",
)


text:
const _PROC_REG_RKEY = "proc_reg"

const _PROC_REG_EXT = ".proc.toml"
_proc_reg_file(gw::GitWorker, pid, rid) = joinpath(gw_procs_dir(gw), string(pid, ".", rid, _PROC_REG_EXT))
_is_proc_reg_file(fn) = endswith(fn, _PROC_REG_EXT)

function _reg_proc!(gw::GitWorker, pid, rid)

    # reg file
    rfile = _proc_reg_file(gw, rid, pid)

    # get proc info
    lstart = _get_proc_lstart(pid)
    isempty(lstart) && return false

    # up ram
    sdat = reg_sdat!(gw, rfile;
        rid, lstart, pid
    )

    # up disk
    _mkdir(rfile)
    _write_toml(rfile, sdat)

    return rfile
end

function _find_proc_reg(gw::GitWorker, khint)
    khint = basename(string(khint))
    pdir = gw_procs_dir(gw)
    _find_by_khint(gw, pdir, khint)
end

function _up_procs_from_disk(gw::GitWorker)
    pdir = gw_procs_dir(gw)
    empty_sreg!(gw, pdir)
    for path in _readdir(pdir; join = true)
        !_is_proc_reg_file(path) && continue
        fdat = _read_toml(path)
        !haskey(fdat, "pid") && continue
        reg_sdat!(gw, path, fdat)
    end
    return gw
end

function _validate_proc(gw::GitWorker, phint)
    pdat = _find_proc_reg(gw, phint)
    isempty(pdat) && return false
    rlstart = get(pdat, "lstart", "")
    isempty(rlstart) && return false
    pid = get(pdat, "pid", "")
    isempty(pid) && return false
    clstart = _get_proc_lstart(pid)
    isempty(clstart) && return false
    return rlstart == clstart
end

max_stars_repo_name: josePereiro/GitWorkers

max_issues_repo_name: josePereiro/GitWorkers

max_forks_repo_name: josePereiro/GitWorkers