text:
#using PyPlot
#pygui(true)
using Statistics
using Random
using Plots
import Gadfly
using DataFrames

mutable struct Agent
    opinion :: AbstractFloat

    function Agent()
        new(rand()*2 -1 )
    end
end

function talkto!(target::Agent, source::Agent, bounded_rat)
    if abs(target.opinion - source.opinion) < bounded_rat
        target.opinion = (target.opinion + source.opinion)/2
    end
end

function simulation(agent_count, steps, bounded_rat = 0.2)
    agents = [Agent() for x in 1:agent_count]

    PyPlot.clf()
    hist([a.opinion for a in agents])
    for step = 1:steps

        @inbounds source = agents[rand(1:agent_count)]
        @inbounds target = agents[rand(1:agent_count)]
        talkto!(target, source, bounded_rat)
    end

    hist([a.opinion for a in agents])
    #display(gcf())
    df = DataFrame(opinions = [a.opinion for a in agents])
    Gadfly.plot(df, x=:opinions, Gadfly.Geom.point)

    (mean([a.opinion for a in agents]), std([a.opinion for a in agents]))
end

@time res = simulation(100000, 10000, 0.4)
println(res)

ext: jl


text:
include("../utility.jl")
using Test
using Distributions, Turing

@model gdemo(x) = begin
  s ~ InverseGamma(2, 3)
  m ~ Normal(0, sqrt(s))
  x[1] ~ Normal(m, sqrt(s))
  x[2] ~ Normal(m, sqrt(s))
  return s, m
end

model_f = gdemo([1.5, 2.0])

alg = NUTS(5000, 1000, 0.65)
res = sample(model_f, alg)

@info(mean(res[:s][1000:end])," ≈ ", 49/24, "?")
@info(mean(res[:m][1000:end])," ≈ ", 7/6, "?")
@test mean(res[:s][1000:end]) ≈ 49/24 atol=0.6
@test mean(res[:m][1000:end]) ≈ 7/6 atol=0.6


text:
# nextpow2(W) = vshl(one(W), vsub_fast(8sizeof(W), leading_zeros(vsub_fast(W, one(W)))))

# @inline _pick_vector(::StaticInt{W}, ::Type{T}) where {W,T} = Vec{W,T}
# @inline pick_vector(::Type{T}) where {T} = _pick_vector(pick_vector_width(T), T)
# @inline function pick_vector(::Val{N}, ::Type{T}) where {N, T}
#     _pick_vector(smin(nextpow2(StaticInt{N}()), pick_vector_width(T)), T)
# end
# pick_vector(N::Int, ::Type{T}) where {T} = pick_vector(Val(N), T)

@inline MM(::Union{Val{W},StaticInt{W}}) where {W} = MM{W}(0)
@inline MM(::Union{Val{W},StaticInt{W}}, i) where {W} = MM{W}(i)
@inline MM(::Union{Val{W},StaticInt{W}}, i::AbstractSIMDVector{W}) where {W} = i
@inline MM(::StaticInt{W}, i, ::StaticInt{X}) where {W,X} = MM{W,X}(i)
@inline gep(ptr::Ptr, i::MM) = gep(ptr, data(i))

@inline Base.one(::Type{MM{W,X,I}}) where {W,X,I} = one(I)
@inline staticm1(i::MM{W,X,I}) where {W,X,I} = MM{W,X}(vsub_fast(data(i), one(I)))
@inline staticp1(i::MM{W,X,I}) where {W,X,I} = MM{W,X}(vadd_fast(data(i), one(I)))
@inline vadd_fast(i::MM{W,X}, j::Integer) where {W,X} = MM{W,X}(vadd_fast(data(i), j))
@inline vadd_fast(i::Integer, j::MM{W,X}) where {W,X} = MM{W,X}(vadd_fast(i, data(j)))
@inline vadd_fast(i::MM{W,X}, ::StaticInt{j}) where {W,X,j} = MM{W,X}(vadd_fast(data(i), StaticInt{j}()))
@inline vadd_fast(::StaticInt{i}, j::MM{W,X}) where {W,X,i} = MM{W,X}(vadd_fast(StaticInt{i}(), data(j)))
@inline vadd_fast(i::MM{W,X}, ::StaticInt{0}) where {W,X} = i
@inline vadd_fast(::StaticInt{0}, j::MM{W,X}) where {W,X} = j
@inline vsub_fast(i::MM{W,X}, j::Integer) where {W,X} = MM{W,X}(vsub_fast(data(i), j))
@inline vsub_fast(i::MM{W,X}, ::StaticInt{j}) where {W,X,j} = MM{W,X}(vsub_fast(data(i), StaticInt{j}()))
@inline vsub_fast(i::MM{W,X}, ::StaticInt{0}) where {W,X} = i

@inline vadd(i::MM{W,X}, j::Integer) where {W,X} = MM{W,X}(vadd_fast(data(i), j))
@inline vadd(i::Integer, j::MM{W,X}) where {W,X} = MM{W,X}(vadd_fast(i, data(j)))
@inline vadd(i::MM{W,X}, ::StaticInt{j}) where {W,X,j} = MM{W,X}(vadd_fast(data(i), StaticInt{j}()))
@inline vadd(::StaticInt{i}, j::MM{W,X}) where {W,X,i} = MM{W,X}(vadd_fast(StaticInt{i}(), data(j)))
@inline vsub(i::MM{W,X}, j::Integer) where {W,X} = MM{W,X}(vsub_fast(data(i), j))
@inline vsub(i::MM{W,X}, ::StaticInt{j}) where {W,X,j} = MM{W,X}(vsub_fast(data(i), StaticInt{j}()))
@inline vsub(i::MM{W,X}, ::StaticInt{0}) where {W,X} = i
@inline vsub(i::MM) = i * StaticInt{-1}()
@inline vsub_fast(i::MM) = i * StaticInt{-1}()
@inline vmul(::StaticInt{M}, i::MM{W,X}) where {M,W,X} = MM{W}(vmul_fast(data(i), StaticInt{M}()), StaticInt{X}() * StaticInt{M}())
@inline vmul(i::MM{W,X}, ::StaticInt{M}) where {M,W,X} = MM{W}(vmul_fast(data(i), StaticInt{M}()), StaticInt{X}() * StaticInt{M}())
@inline vrem(i::MM{W,X,I}, ::Type{I}) where {W,X,I<:IntegerTypesHW} = i
@inline vrem(i::MM{W,X}, ::Type{I}) where {W,X,I<:IntegerTypesHW} = MM{W,X}(data(i) % I)
@inline veq(::AbstractIrrational, ::MM{W,<:Integer}) where {W} = zero(Mask{W})
@inline veq(x::AbstractIrrational, i::MM{W}) where {W} = x == Vec(i)
@inline veq(::MM{W,<:Integer}, ::AbstractIrrational) where {W} = zero(Mask{W})
@inline veq(i::MM{W}, x::AbstractIrrational) where {W} = Vec(i) == x

@inline function vsub(i::NativeTypes, j::MM{W,X}) where {W,X}
    MM(StaticInt{W}(), vsub(i, data(j)), -StaticInt{X}())
end
@inline function vsub_fast(i::NativeTypes, j::MM{W,X}) where {W,X}
    MM(StaticInt{W}(), vsub_fast(i, data(j)), -StaticInt{X}())
end

@inline function Base.in(m::MM{W,X,<:Integer}, r::AbstractUnitRange) where {W,X}
    vm = Vec(m)
    (vm ≥ first(r)) & (vm ≤ last(r))
end

@inline function pick_vector_width_shift(args::Vararg{Any,K}) where {K}
    W = pick_vector_width(args...)
    W, intlog2(W)
end


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

using Markdown
using InteractiveUtils

# ╔═╡ 22a37ff4-3ae0-4338-86a2-3310b7773c1d
begin
    import Pkg
    # activate the shared project environment|
    Pkg.activate(Base.current_project())
    using Omega
	using Distributions
	using UnicodePlots
	using Random: MersenneTwister
end

# ╔═╡ 6b2fe98e-76cf-11eb-39f8-9ff52d23e688
md"# Bayesian Linear Regression"

# ╔═╡ a2fb35f7-8c48-4e70-aff8-d8ee431aa872
md"Bayesian linear regression is a method to infer the parameters of linear model between two variables.  The \"Bayesian\" part of Bayesian linear regression means means that it treats the problem of finding these parameters as one of Bayesian inference."

# ╔═╡ b583b52a-76ce-11eb-088b-dbeb7791729f
md"""## Data
First, let's create some fake data"""

# ╔═╡ ee2f8174-76ce-11eb-2c74-298d6072afee
rng = MersenneTwister(0);

# ╔═╡ d10a9c5a-76ce-11eb-0ece-bd1b6e946fe0
N = 100;

# ╔═╡ d83b52d0-76ce-11eb-0d3f-d1ddd64d6e08
xs = rand(rng, N)

# ╔═╡ e61b99a8-76ce-11eb-2d0e-1f8d2f31eba4
m = 0.8

# ╔═╡ 0af4f988-76cf-11eb-012e-cf2182fa826b
c = 1.2

# ╔═╡ ff7d2e7c-76ce-11eb-249b-cbf84db4fe9b
linear_model(x, m, c) = m * x + c

# ╔═╡ 2288d0d2-76d0-11eb-29e9-21a78bc370e8
obs_model(x) = linear_model(x, m, c) + randn(rng) * 0.1;

# ╔═╡ 107a8c36-76cf-11eb-1adf-43f11a9a4204
ys = obs_model.(xs);

# ╔═╡ 1dc0257e-76cf-11eb-3fa5-377113ba0a49
UnicodePlots.scatterplot(xs, ys)

# ╔═╡ 61f89028-76cf-11eb-1813-9179e3467f47
md"## Probabilistic Model"

# ╔═╡ 79cd6e26-76cf-11eb-1930-399093becba8
M = @~ Normal(0, 1);

# ╔═╡ a9765e1c-76cf-11eb-3495-79aa95e21afd
C = @~ Normal(0, 1);

# ╔═╡ 7fce6a74-8b68-45ce-bb9c-a4c4e39fcc1e
struct ϵ end

# ╔═╡ b2cc4db2-76cf-11eb-390d-f151ff4b939d
Y_class(i, ω) = linear_model(xs[i], M(ω), C(ω)) + (ϵ ∘ i ~ Normal(0, 0.1))(ω);

# ╔═╡ 93ed05dd-cc5a-4b77-af0d-7cca48777510
revise

# ╔═╡ 800049e2-f44b-4839-a2ae-7fb0beeb5c0e
Y⃗ = Mv(1:N, Y_class)

# ╔═╡ fa8b5cc2-7728-11eb-24b3-5db1bd7d0255
UnicodePlots.scatterplot(xs, randsample(Y⃗))

# ╔═╡ cf09331c-c5cd-49f0-8ad5-eced9cb56d37
(ω = defω(); Y⃗(ω); ω)

# ╔═╡ b3c0e452-fff8-419e-945c-27e2d21c9fe3
evidence = pw(==, Y⃗, ys)

# ╔═╡ 3d5b3ac2-bf70-4710-9c97-1cbef7be4c11
joint_posterior = @joint(M, C) |ᶜ evidence

# ╔═╡ aa781fbe-76d0-11eb-387a-3dc1cbf700f8
samples = randsample(joint_posterior |ᶜ evidence, 1000; alg = MH)

# ╔═╡ f7efb0d0-76d1-11eb-3440-47df94aeb74e
UnicodePlots.scatterplot(xs, samples)

# ╔═╡ Cell order:
# ╟─6b2fe98e-76cf-11eb-39f8-9ff52d23e688
# ╟─a2fb35f7-8c48-4e70-aff8-d8ee431aa872
# ╠═22a37ff4-3ae0-4338-86a2-3310b7773c1d
# ╟─b583b52a-76ce-11eb-088b-dbeb7791729f
# ╠═ee2f8174-76ce-11eb-2c74-298d6072afee
# ╠═d10a9c5a-76ce-11eb-0ece-bd1b6e946fe0
# ╠═d83b52d0-76ce-11eb-0d3f-d1ddd64d6e08
# ╠═e61b99a8-76ce-11eb-2d0e-1f8d2f31eba4
# ╠═0af4f988-76cf-11eb-012e-cf2182fa826b
# ╠═ff7d2e7c-76ce-11eb-249b-cbf84db4fe9b
# ╠═2288d0d2-76d0-11eb-29e9-21a78bc370e8
# ╠═107a8c36-76cf-11eb-1adf-43f11a9a4204
# ╠═1dc0257e-76cf-11eb-3fa5-377113ba0a49
# ╟─61f89028-76cf-11eb-1813-9179e3467f47
# ╠═79cd6e26-76cf-11eb-1930-399093becba8
# ╠═a9765e1c-76cf-11eb-3495-79aa95e21afd
# ╠═7fce6a74-8b68-45ce-bb9c-a4c4e39fcc1e
# ╠═b2cc4db2-76cf-11eb-390d-f151ff4b939d
# ╠═93ed05dd-cc5a-4b77-af0d-7cca48777510
# ╠═800049e2-f44b-4839-a2ae-7fb0beeb5c0e
# ╠═fa8b5cc2-7728-11eb-24b3-5db1bd7d0255
# ╠═cf09331c-c5cd-49f0-8ad5-eced9cb56d37
# ╠═b3c0e452-fff8-419e-945c-27e2d21c9fe3
# ╠═3d5b3ac2-bf70-4710-9c97-1cbef7be4c11
# ╠═aa781fbe-76d0-11eb-387a-3dc1cbf700f8
# ╠═f7efb0d0-76d1-11eb-3440-47df94aeb74e

max_stars_repo_name: zenna/expect

max_issues_repo_name: zenna/expect

max_forks_repo_name: zenna/expect


text:
# add one unit to the voltage setting
# SV+
export psp_incrementvoltage

const INCREMENT_VOLTAGE = [0x53, 0x56, 0x2B, 0x0D]

"""
    psp_incrementvoltage(io)

Add one unit to the voltage setting.
"""
function psp_incrementvoltage(io_psp::IO)
  write(io_psp, INCREMENT_VOLTAGE)
  return nothing
end


text:
module PH2514

import Base
# Write your package code here.
export C, λ_cf, f_cλ, f_hE, E_hf, E_Zn, Energy, ElectronVolt, Joule, eV, J, Wavelength, Frequency, WiensDisplacementLaw, SteffannBoltzman, Vrms, GyroRadius, GyroFrequency, PlasmaFrequency

C = Dict(
    :c => 299792458,
    :h => 6.626e-34,
    :J2eV => 6.242e18,
    :eV2J => 1.602e-19,
    :E1 => -13.58,
    :k => 1.380649e-23,
    :σ => 5.67e-8,
    :a =>2.898e-3,
    :me => 9.10938356e-31,
    :mp => 1.6726e-27,
    :qe => 1.60217662e-19,
    :ϵ0 => 8.85418e-12,
)

abstract type Unit <: Number end
abstract type Energy <: Unit end
struct ElectronVolt <: Energy
    value::Number
end
struct Joule <: Energy
    value::Number
end
struct Wavelength <: Unit
    value::Number
end
struct Frequency <: Unit
    value::Number
end
promote_rule(::Type{Number},    ::Type{Unit})           = Unit
promote_rule(::Type{Unit},      ::Type{Number})         = Unit
promote_rule(::Type{Joule},     ::Type{ElectronVolt})   = Joule
promote_rule(::Type{Number},    ::Type{Frequency})      = Frequency
promote_rule(::Type{Frequency}, ::Type{Number})         = Frequency
promote_rule(::Type{Number},    ::Type{Wavelength})     = Wavelength
promote_rule(::Type{Wavelength},::Type{Number})         = Wavelength

ElectronVolt(J::Joule)  = convert(ElectronVolt,J)
ElectronVolt(f::Frequency) = convert(ElectronVolt,f)
ElectronVolt(λ::Wavelength) = convert(ElectronVolt,λ)
Joule(eV::ElectronVolt) = convert(Joule,eV)
Joule(f::Frequency) = convert(Joule,f)
Joule(λ::Wavelength) = convert(Joule,λ)

Wavelength(f::Frequency) = convert(Wavelength,f)
Frequency(λ::Wavelength) = convert(Frequency,λ)
Wavelength(e::Energy)    = convert(Wavelength,e)
Frequency(e::Energy)     = convert(Frequency,e)

Base.convert(::Type{ElectronVolt},  t::Joule)           = ElectronVolt(t.value * C[:J2eV])
Base.convert(::Type{Joule},         t::ElectronVolt)    = Joule(t.value * C[:eV2J])

Base.convert(::Type{Wavelength},    t::Frequency)       = Wavelength(C[:c] / t.value)
Base.convert(::Type{Frequency},     t::Wavelength)      = Frequency(C[:c] / t.value)

Base.convert(::Type{Frequency},     t::Energy)          = Frequency(convert(Joule,t).value / C[:h])
Base.convert(::Type{Joule},         t::Frequency)       = Joule(C[:h] * t.value)
Base.convert(::Type{ElectronVolt},  t::Frequency)       = ElectronVolt(convert(Joule,t))
Base.convert(::Type{Joule},         t::Wavelength)      = convert(Joule,Frequency(t))
Base.convert(::Type{ElectronVolt},  t::Wavelength)      = convert(ElectronVolt,Frequency(t))

Base.convert(::Type{Number},        t::Unit)            = Number(t.value)
Base.convert(::Type{Unit},          t::Number)          = Unit(t)
Base.convert(::Type{Wavelength},    t::Number)          = Wavelength(t)
Base.convert(::Type{Frequency},     t::Number)          = Frequency(t)

Base. +(a::Joule,b::Joule) = Joule(a.value + b.value)
Base. -(a::Joule,b::Joule) = Joule(a.value + b.value)
Base. +(a::ElectronVolt,b::ElectronVolt) = ElectronVolt(a.value + b.value)
Base. -(a::ElectronVolt,b::ElectronVolt) = ElectronVolt(a.value + b.value)
Base. +(a::Wavelength,b::Wavelength) = Wavelength(a.value + b.value)
Base. -(a::Wavelength,b::Wavelength) = Wavelength(a.value + b.value)
Base. +(a::Frequency,b::Frequency) = Frequency(a.value + b.value)
Base. -(a::Frequency,b::Frequency) = Frequency(a.value + b.value)

Base. -(a::Unit,b::Number) = -(promote(a,b)...)
Base. +(a::Unit,b::Number) = +(promote(a,b)...)
Base. *(a::Unit,b::Number) = *(promote(a,b)...)
Base. ÷(a::Unit,b::Number) = ÷(promote(a,b)...)
Base. *(a::Unit,b::Unit)   = *(promote(a,b)...)

function Base.show(io::IO, J::Joule)
    print(io, "$(J.value) J")
end

function Base.show(io::IO, e::ElectronVolt)
    print(io, "$(e.value) eV")
end

function Base.show(io::IO, λ::Wavelength)
    print(io, "$(λ.value) m")
end

function Base.show(io::IO, f::Frequency)
    print(io, "$(f.value) Hz")
end

const J = Joule
const eV = ElectronVolt

λ_cf(f::Frequency) = Wavelength(f)

f_cλ(λ::Wavelength) = Frequency(λ)

f_hE(E::Energy) = Frequency(E)

E_hf(f::Frequency) = Joule(f)

function E_Zn(Z,n)
    C[:E1] * (Z/n)^2
end

function WiensDisplacementLaw(T)
    C[:a]/T
end

function PlanckRadiation_L(λ,T)
    2 * C[:h] * C[:c]^2 / λ^5 / ( exp(C[:h]*C[:c]/(λ*C[:k]*T)) - 1)
end

function SteffannBoltzman(ϵ,T)
    C[:σ] * ϵ * T^4
end

function Vrms(m,T)
    √(3 * C[:k] * T / m)
end

function GyroRadius(m,v,q,B)
    (m*v) / (q*B)
end

function GyroFrequency(q,B,m)
    q*B/m
end

function PlasmaFrequency(n)
    √((n * C[:qe]^2)/(C[:ϵ0]*C[:me]))
end

end


text:
using Flux.Optimise
using Flux.Optimise: runall
using Flux: Params, gradient
using Test

@testset "Optimise" begin
  w = randn(10, 10)
  @testset for opt in [ADAMW(), ADAGrad(0.1), AdaMax(), ADADelta(0.9), AMSGrad(),
                       NADAM(), RADAM(), Descent(0.1), ADAM(), Nesterov(), RMSProp(),
                       Momentum()]
    w′ = randn(10, 10)
    loss(x) = Flux.mse(w*x, w′*x)
    for t = 1: 10^5
      θ = Params([w′])
      x = rand(10)
      θ̄ = gradient(() -> loss(x), θ)
      Optimise.update!(opt, θ, θ̄)
    end
    @test loss(rand(10, 10)) < 0.01
  end
end

@testset "Optimiser" begin
  w = randn(10, 10)
  @testset for Opt in [InvDecay, WeightDecay, ExpDecay]
    w′ = randn(10, 10)
    loss(x) = Flux.mse(w*x, w′*x)
    opt = Optimiser(Opt(), ADAM(0.001))
    for t = 1:10^5
      θ = Params([w′])
      x = rand(10)
      θ̄ = gradient(() -> loss(x), θ)
      Optimise.update!(opt, θ, θ̄)
    end
    @test loss(rand(10, 10)) < 0.01
  end
end

@testset "Training Loop" begin
  i = 0
  l = 1

  Flux.train!(() -> (sleep(0.1); i += 1; l),
              (),
              Iterators.repeated((), 100),
              Descent(),
              cb = Flux.throttle(() -> (i > 3 && Flux.stop()), 1))

  @test 3 < i < 50

  # Test multiple callbacks
  x = 0
  fs = [() -> (), () -> x = 1]
  cbs = runall(fs)
  cbs()
  @test x == 1
end

@testset "ExpDecay" begin
    w = randn(10, 10)
    o = ExpDecay(0.1, 0.1, 1000, 1e-4)
    w1 = randn(10,10)
    loss(x) = Flux.mse(w*x, w1*x)
    flag = 1
    decay_steps = []
    for t = 1:10^5
      prev_eta = o.eta
      θ = Params([w1])
      x = rand(10)
      θ̄ = gradient(() -> loss(x), θ)
      prev_grad = collect(θ̄[w1])
      delta = Optimise.apply!(o, w1, θ̄[w1])
      w1 .-= delta
      new_eta = o.eta
      if new_eta != prev_eta
        push!(decay_steps, t)
      end
      array = fill(o.eta, size(prev_grad))
      if array .* prev_grad != delta
        flag = 0
      end
    end
    @test flag == 1
    # Test to check if decay happens at decay steps. Eta reaches clip value eventually.
    ground_truth = []
    for i in 1:11
      push!(ground_truth, 1000*i)  # Expected decay steps for this example.
    end
    @test decay_steps == ground_truth
    @test o.eta == o.clip
end


text:
cd(@__DIR__); import Pkg; Pkg.activate(".")
coords = (left = -145.0, right = -50.0, bottom = 20.0, top = 75.0)
Pkg.develop(path="$(homedir())/github/SimpleSDMLayers.jl")
using Revise
using SimpleSDMLayers
using Test

## Step 8 : DataFrames no nothing ####

# For issue
using SimpleSDMLayers
using DataFrames

layer = SimpleSDMPredictor(WorldClim, BioClim, 1)

df = DataFrame([layer, layer])
allowmissing!(df)
for col in [:x1, :x2]
    replace!(df[!, col], nothing => missing)
end
dropmissing(df, [:x1, :x2])

# Tests
DataFrame(layer)
DataFrame([layer, layer])

l1 = layer[top = -89.0]
typeof(l1)
eltype(l1)
typeof(l1.grid)
eltype(l1.grid)

df1 = DataFrame(l1)
typeof(df1.values)
eltype(df1.values)

layers = [layer, layer]
layers = [temperature_clip, temperature_clip]


text:
using RankAggregation
using DataFrames
using Test

@testset "RankAggregation.jl" begin
  objects = DataFrame(object=[:a,:b,:c],
                      score1=[0.9, 0.7, 0.5],
                      score2=[0.8, 0.9, 0.4])
  @test rank(objects, :score1) == [1,2,3]
  @test rank(objects, :score2) == [2,1,3]
  @test rank(objects, (:score1,:score2)) == [1,2,3]
  @test rank(objects, :score2, rev=true) == [2,3,1]
end


text:
import StaticArrays: SVector

@testset "Newton" begin
    @polyvar x y z

    f = x^2*y + y^2 + 3x^2*y^3 + y^4 + x^4*y^4

    AF = AmoebaFiber2D(f, (0, 0))
    PolynomialAmoebas.update_fiber!(AF, (0, 1))
    @test PolynomialAmoebas.newton(AF, SVector(0.5π, 1.5π), 100, 1e-6).converged

    # underdetermined
    AF = AmoebaFiber3D(x^2+y^2+z^2+1, (0, 1, 0))
    PolynomialAmoebas.update_fiber!(AF, (0, 0, 0))
    v = SVector(3.97394, 2.84526, 1.69766)
    @test PolynomialAmoebas.newton(AF, v, 150, 1e-7).converged
end


text:
module Static

export serve

using LiveServer
using ..Types
using ..Parse
using ..Write
using ..BuildTools
using ..BuildTools: skipfile

function serve(p::String=pwd(); kwargs...)
    serve(parse_slide(p); kwargs...)
end

function serve(p::Presentation; kwargs...)
    build_dir = build(p)
    raw_files = collect_files(p)
    assets = assets_files(p)
    # TODO: update file sets
    function reload_slide(fp)
        # TODO: only build changed file
        if fp in assets
            @info "asset $fp changed, update asset"
            cp(fp, joinpath(build_dir, fp); force=true)
        end

        if fp in raw_files
            @info "content $fp changed, recompile index.html"
            write(joinpath(build_dir, "index.html"), parse_slide(p.path))
        end

        LiveServer.file_changed_callback(fp)
    end

    fw = LiveServer.SimpleWatcher(reload_slide)

    for file in raw_files
        LiveServer.watch_file!(fw, file)
    end

    for file in assets
        LiveServer.watch_file!(fw, file)
    end

    return LiveServer.serve(fw; dir=build_dir, kwargs...)
end

function collect_files(p::Presentation)
    if p.standalone
        return [p.path]
    else
        files = [p.path]
        for slide in p.slides
            append!(files, collect_files(slide))
        end
        return unique(files)
    end
end

collect_files(p::BasicSlide) = [p.path]

function collect_files(p::SectionSlide)
    files = [p.path]
    for slide in p.slides
        append!(files, collect_files(slide))
    end

    return unique(files)
end

function assets_files(p::Presentation)
    p.standalone && return []

    files = []
    for dir in readdir(p.path)
        if !startswith(dir, "_") && !startswith(dir, ".") && !skipfile(dir)
            if isfile(dir)
                push!(files, dir)
            else
                for (root, dirs, files) in walkdir(dir)
                    for file in files
                        push!(files, joinpath(root, file))
                    end
                end
            end
        end
    end
    return files
end

end


text:
# Aggregate functions.

mutable struct AggregateClause <: AbstractSQLClause
    name::Symbol
    distinct::Bool
    args::Vector{SQLClause}
    filter::Union{SQLClause, Nothing}
    over::Union{SQLClause, Nothing}

    AggregateClause(;
                    name::Union{Symbol, AbstractString},
                    distinct = false,
                    args = SQLClause[],
                    filter = nothing,
                    over = nothing) =
        new(Symbol(name), distinct, args, filter, over)
end

AggregateClause(name; distinct = false, args = SQLClause[], filter = nothing, over = nothing) =
    AggregateClause(name = name, distinct = distinct, args = args, filter = filter, over = over)

AggregateClause(name, args...; distinct = false, filter = nothing, over = nothing) =
    AggregateClause(name, distinct = distinct, args = SQLClause[args...], filter = filter, over = over)

"""
    AGG(; name, distinct = false, args = [], filter = nothing, over = nothing)
    AGG(name; distinct = false, args = [], filter = nothing, over = nothing)
    AGG(name, args...; distinct = false, filter = nothing, over = nothing)

An application of an aggregate function.

# Examples

```jldoctest
julia> c = AGG(:COUNT, OP("*"));

julia> print(render(c))
COUNT(*)
```

```jldoctest
julia> c = AGG(:COUNT, distinct = true, :year_of_birth);

julia> print(render(c))
COUNT(DISTINCT "year_of_birth")
```

```jldoctest
julia> c = AGG(:COUNT, OP("*"), filter = OP(">", :year_of_birth, 1970));

julia> print(render(c))
(COUNT(*) FILTER (WHERE ("year_of_birth" > 1970)))
```

```jldoctest
julia> c = AGG(:ROW_NUMBER, over = PARTITION(:year_of_birth));

julia> print(render(c))
(ROW_NUMBER() OVER (PARTITION BY "year_of_birth"))
```
"""
AGG(args...; kws...) =
    AggregateClause(args...; kws...) |> SQLClause

dissect(scr::Symbol, ::typeof(AGG), pats::Vector{Any}) =
    dissect(scr, AggregateClause, pats)

function PrettyPrinting.quoteof(c::AggregateClause, ctx::QuoteContext)
    ex = Expr(:call, nameof(AGG), string(c.name))
    if c.distinct
        push!(ex.args, Expr(:kw, :distinct, c.distinct))
    end
    append!(ex.args, quoteof(c.args, ctx))
    if c.filter !== nothing
        push!(ex.args, Expr(:kw, :filter, quoteof(c.filter, ctx)))
    end
    if c.over !== nothing
        push!(ex.args, Expr(:kw, :over, quoteof(c.over, ctx)))
    end
    ex
end

rebase(c::AggregateClause, c′) =
    AggregateClause(name = c.name, distinct = c.distinct, args = c.args, filter = c.filter, over = rebase(c.over, c′))


text:
function _make_contact_gui()

    grid=Grid()

    training_frame = Frame("Training Data")
    training_grid = Grid()
    push!(training_frame,training_grid)
    grid[1,1]=training_frame

    training_grid[1,1]=Label("Total Training Frames: ")
    training_num_label=Label("1")
    training_grid[2,1]=training_num_label

    tracked_frame = Frame("Predicted Data")
    tracked_grid = Grid()
    push!(tracked_frame,tracked_grid)
    grid[1,2]=tracked_frame

    tracked_load_button = Button("Load Predicted Contact and Pole")
    tracked_grid[1,1] = tracked_load_button

    classifier_frame = Frame("Classifier")
    classifier_grid = Grid()
    push!(classifier_frame, classifier_grid)
    grid[1,3]=classifier_frame

    predictors_box=Box(:v)
    push!(predictors_box,Label("Predictors: "))

    pred_pole_button = CheckButton("Distance to Pole")
    push!(predictors_box,pred_pole_button)

    pred_pole_position = CheckButton("Pole Position")
    push!(predictors_box,pred_pole_position)

    pred_curv = CheckButton("Curvature")
    push!(predictors_box,pred_curv)

    classifier_grid[1,1]=predictors_box

    fit_button = Button("Fit Classifier")
    classifier_grid[1,2]=fit_button

    n_estimators_button = SpinButton(1:1000)
    setproperty!(n_estimators_button,:value,100)
    classifier_grid[2,2] = n_estimators_button

    classifier_grid[3,2] = Label("Number of Estimators")

    forest_depth_button = SpinButton(1:100)
    setproperty!(forest_depth_button,:value,10)
    classifier_grid[4,2] = forest_depth_button

    classifier_grid[5,2] = Label("Forest Depth")

    classifier_grid[1,3] = Label("Cross-validation prediction accuracy: ")
    cv_label = Label("")
    classifier_grid[2,3] = cv_label

    win = Window(grid)
    Gtk.showall(win)
    visible(win,false)

    c_widgets=contact_widgets(win,training_num_label,fit_button,tracked_load_button,
    n_estimators_button,forest_depth_button,cv_label,pred_pole_button,
    pred_pole_position,pred_curv)
end

function add_contact_callbacks(w,handles)

    signal_connect(contact_fit_cb,w.fit_button,"clicked",Void,(),false,(handles,))
    signal_connect(contact_load_predicted_cb,w.load_predicted_button,"clicked",Void,(),false,(handles,))

    nothing
end

function contact_fit_cb(w::Ptr,user_data::Tuple{Tracker_Handles})

    han, = user_data

    han.class.clf=RandomForestClassifier(n_estimators=han.class.n_estimators,max_depth=han.class.forest_depth,random_state=0)

    #Select which variables to use
    if getproperty(han.c_widgets.pred_pole_button,:active,Bool)

    end

    if getproperty(han.c_widgets.pred_pole_position,:active,Bool)

    end

    if getproperty(han.c_widgets.pred_curv,:active,Bool)

    end

    #Remove NaNs
    #han.class.predictors[isnan.(han.class.predictors)] .= 0.0

    #Create prediction
    #ScikitLearn.fit!(han.class.clf,han.class.predictors,han.tracked_contact)

    #Cross validation
    #cross_val_rf=mean(cross_val_score(han.class.clf, han.class.predictors, han.tracked_contact, cv=10))

    #setproperty!(han.c_widgets.cv_label,:label,han.class.cv)

    nothing
end

#Cross validate function

#Predict Function

contact_estimators(han)=getproperty(han.c_widgets.n_estimators_button,:value,Int64)
contact_depth(han)=getproperty(han.c_widgets.forest_depth_button,:value,Int64)

#Load labels
function load_contact_cb(w::Ptr,user_data::Tuple{Tracker_Handles})

    han, = user_data

    filepath = open_dialog("Load Contact Data",han.b["win"])

    if filepath != ""
        file = matopen(filepath,"r")
        if MAT.exists(file,"Touch_Inds")
            touch_inds = read(file,"Touch_Inds")
            han.touch_frames_i = touch_inds
            setproperty!(han.contact_widgets.training_num_label,:label,string(length(han.touch_frames_i)))
        end
        if MAT.exists(file,"Touch")
            touch = read(file,"Touch")
            han.touch_frames = convert(BitArray{1},touch)
        end
        close(file)
    end

    nothing
end

function save_contact_cb(w::Ptr,user_data::Tuple{Tracker_Handles})

    han, = user_data

    filepath = save_dialog("Save Contact Data",han.b["win"])

    if filepath != ""
        if filepath[end-3:end]==".mat"
        else
            filepath=string(filepath,".mat")
        end

        file=matopen(filepath,"w")
        write(file,"Touch_Inds",han.touch_frames_i)
        write(file,"Touch",convert(Array{Int64,1},han.touch_frames))
        close(file)
    end

    nothing
end

function contact_load_predicted_cb(w::Ptr,user_data::Tuple{Tracker_Handles})

    han, = user_data

    #Load Contact
    filepath = open_dialog("Load Predicted Labels",han.b["win"])

    if filepath != ""

        file = matopen(filepath,"r")
        if MAT.exists(file,"Contact")
            contact = read(file,"Contact")
            han.tracked_contact = contact
        end
    end

    filepath = open_dialog("Load Tracked Pole",han.b["win"])

    if filepath != ""

        p=read_pole_hdf5(filepath)
        han.tracked_pole = p
    end

    nothing
end

# Mark Contact Training Frame or not training frame
function touch_override_cb(w::Ptr,user_data::Tuple{Tracker_Handles,Int64})

    han, contact = user_data

    #Check if this frame is already in the contact list
    if isempty(find(han.touch_frames_i.==han.displayed_frame))

        push!(han.touch_frames_i,han.displayed_frame)
        push!(han.touch_frames,contact)

        setproperty!(han.contact_widgets.training_num_label,:label,string(length(han.touch_frames_i)))

        #Do we need to sort these here?

    else #Change existing entry

        ind=find(han.touch_frames_i.==han.displayed_frame)[1]
        han.touch_frames[ind] = contact
    end

    draw_touch(han)

    nothing
end

#=
Draw marker to indicate that touch has occured
=#

function draw_touch(han::Tracker_Handles)

    if !isempty(find(han.touch_frames_i.==han.displayed_frame))

        ctx=Gtk.getgc(han.c)

        ind=find(han.touch_frames_i.==han.displayed_frame)[1]

        if han.touch_frames[ind]
            set_source_rgb(ctx,1,0,0)
        else
            set_source_rgb(ctx,1,1,1)
        end

        rectangle(ctx,600,0,20,20)
        fill(ctx)

        reveal(han.c)
    end

    nothing
end

function draw_touch_prediction(han::Tracker_Handles)

    if (han.show_contact)

        try
            ctx = Gtk.getgc(han.c)

            w = 640
            h = 480

            if han.tracked_contact[han.displayed_frame]
                set_source_rgb(ctx,1,0,0)
            else
                set_source_rgb(ctx,1,1,1)
            end

            rectangle(ctx, w - 40, h - 40, 20, 20)
            fill(ctx)

            reveal(han.c)
        catch
            println("Could not draw contact prediction")
        end
    end

    nothing
end


text:
#=
block:
- Julia version:
- Author: miros
- Date: 2019-01-16
=#

using Distributed

# @everywhere begin
    using SparseArrays
    import SparseArrays.spzeros
    # include("print_function.jl")
    ArrayOrTuple = Union{Array{Int, 1}, Tuple}
    ArrayOrTupleOrNothing = Union{Array{Int, 1}, Tuple, Nothing}

    function get_global_block_size_if_necessary(block_size::ArrayOrTupleOrNothing)
        if block_size == nothing
            return _b3_size
        else
            return block_size
        end
    end

    """
    seg3d
    """
    function number_of_blocks_per_axis(seg3d_size, block_size)
    #     println("block_size: ", block_size)
        dim = nothing
        if typeof(block_size) == Tuple{Int64}
            dim = nfields(seg3d_size)
        elseif  typeof(block_size) == Array{Int64, 1}
            dim = length(block_size)
        else
            tp = typeof(block_size)
            @warn "Unknown type of block_size: $tp"
        end
    #     println( "dim:", dim, " seg size:", seg3d_size, " block size: ", block_size)

        blocks_number = Array{Int}(undef, dim)
        for k in 1:dim
            # print(k)
            delenec = seg3d_size[k]
            delitel = block_size[k]
            number_for_this_axis = cld(delenec, delitel)
            # print(" ", delenec, ", ", delitel, ", ", number_for_this_axis)

            blocks_number[k] = number_for_this_axis
        end
        return prod(blocks_number), blocks_number
    end

    function block_getter_fast(data3d, block_size; fixed_block_size=false)
        if fixed_block_size == false
            @error "block_getter_fast() work just for fixed_block_size"
            exit()
        end
        data_size = size_as_array(size(data3d))
        n, blocks_per_axis = number_of_blocks_per_axis(data_size, block_size)
        return n, (data3d, block_size, 0, blocks_per_axis, fixed_block_size)
    end

    function get_block_fast(block_i::Integer, data3d, block_size:: Array{Int64, 1}, margin_size::Int, blocks_number_axis, fixed_block_size::Bool=false)
    end

    function block_getter(data3d, block_size; fixed_block_size=false)
        data_size = size_as_array(size(data3d))
        n, blocks_per_axis = number_of_blocks_per_axis(data_size, block_size)
        return n, (data3d, block_size, 0, blocks_per_axis, fixed_block_size)
    end

    function get_block(block_i::Integer, data3d, block_size:: Array{Int64, 1}, margin_size::Int, blocks_number_axis, fixed_block_size::Bool=false)
        # This old implementation is slower.
        a = Array{Int}(
            undef,
            blocks_number_axis[1],
            blocks_number_axis[2],
            blocks_number_axis[3]
        )
        bsub = CartesianIndices(a)[block_i]
        bsub_arr = [bsub[1], bsub[2], bsub[3]]

        # TODO use faster fallowing function
        # bsub2 = grid_voxel_id_to_cart(blocks_number_axis, block_i)
        # bsub_arr2 = [bsub2[3], bsub2[2], bsub2[1]]
        # TODO fix the problem
        # if bsub_arr != bsub_arr2
        #     @error "bsub != bsub2", bsub_arr, bsub_arr2, blocks_number_axis, block_i
        #     exit()
        # else
        #     @info "$bsub_arr != $bsub_arr2 " bsub_arr, bsub_arr2
        # end

        data_size = size_as_array(size(data3d))

        first = (bsub_arr .== [1, 1, 1])
        last = (bsub_arr .== blocks_number_axis)
    #     print(" bsub :", bsub, "block number axis", blocks_number_axis, " first last ", first, last, "\n")
        if any(first) || any(last)
            if margin_size == 0
    #             println("running margin_size == 0 version")
                xst, xsp, yst, ysp, zst, zsp = data_sub_from_block_sub(
                    block_size, margin_size, bsub
                )
                stop = [xsp, ysp, zsp]
    #             print("stop ", stop)
                hard_max!(stop, data_size)
            #     new_stop = zeros(eltype(stop), size(stop)...)
    #             print("stop ", stop)
                oxst = xst
                oyst = yst
                ozst = zst
                oxsp, oysp, ozsp = stop
                outdata = data3d[oxst:oxsp, oyst:oysp, ozst:ozsp]

                if fixed_block_size
                    outdata = zeros(
                        eltype(data3d),
                        block_size...
                    )
                    outdata[1:oxsp-oxst+1, 1:oysp-oyst+1, 1:ozsp-ozst+1] = data3d[oxst:oxsp, oyst:oysp, ozst:ozsp]
                else
                    outdata = data3d[oxst:oxsp, oyst:oysp, ozst:ozsp]
                end
            else
    #             print(" end of col, row or slice ", bsub, "\n")

        #         outdata = zeros(
        #             eltype(data3d),
        #             block_size[1] + margin_size,
        #             block_size[2] + margin_size,
        #             block_size[3] + margin_size
        #         )
                xst, xsp, yst, ysp, zst, zsp = data_sub_from_block_sub(
                    block_size, margin_size, bsub
                )
    #             "subs before: "
    #             print_slice3(xst, xsp, yst, ysp, zst, zsp)
                xst, oxst, xsh = get_start_and_outstart_ind(xst, margin_size)
                yst, oyst, ysh = get_start_and_outstart_ind(yst, margin_size)
                zst, ozst, zsh = get_start_and_outstart_ind(zst, margin_size)
        #         print("[", xsh, ", ", ysh, ",", zsh, "]")
                szx, szy, szz = size(data3d)
                bszx, bszy, bszz = block_size
                xsp, oxsp = get_end_and_outend_ind(xst, xsp, szx, xsh)
                ysp, oysp = get_end_and_outend_ind(yst, ysp, szy, ysh)
                zsp, ozsp = get_end_and_outend_ind(zst, zsp, szz, zsh)
    #             print(" sh[",xsh, ",", ysh, ",", zsh, "]")
    #             print("postprocessing input")
                # print_slice3(xst, xsp, yst, ysp, zst, zsp)
    #             print("postprocessing output")
                # print_slice3(oxst, oxsp, oyst, oysp, ozst, ozsp)
        #         print_slice3(xst, xsp, yst, ysp, zst, zsp)
        #         print_slice3(oxst, oxsp, oyst, oysp, ozst, ozsp)
                if fixed_block_size
                    temp_block_size = block_size
                else
                    temp_block_size = [
                        xsh + xsp - xst + 1,
                        ysh + ysp - yst + 1,
                        zsh + zsp - zst + 1
                    ]
                end

                # println("temp_block_size")
                # display(temp_block_size)

                outdata = zeros(
                    eltype(data3d),
                    temp_block_size...
                )
                # println("out data shape")
                # display(size(outdata))
                outdata[oxst:oxsp, oyst:oysp, ozst:ozsp] = data3d[
                    xst:xsp, yst:ysp, zst:zsp
                ]
            end

        else
            # println("not on edge")
            oxst, oxsp, oyst, oysp, ozst, ozsp = data_sub_from_block_sub(
                block_size, margin_size, bsub
            )
    #         print(oxst, ":", oxsp, ", ", oyst, ":", oysp, ", ", ozst, ":", ozsp)
            outdata = data3d[oxst:oxsp, oyst:oysp, ozst:ozsp]
            offset = [oxst, oyst, ozst]
        end
        offset = [oxst - 1, oyst - 1, ozst - 1]
        sz = [oxsp - oxst + 1, oysp - oyst + 1, ozsp - ozst + 1]

        return outdata, offset, sz
    end

    function get_start_and_outstart_ind(xst, margin_size)
        if xst < 1
            oxst = 2 - xst
            xst = 1
            xshift = oxst - xst
        else
            oxst = 1
            xshift = 0
            # xst = xst
        end
        return xst, oxst, xshift
    end

    function get_end_and_outend_ind(xst, xsp, szx, xsh)
        if szx < xsp
    #         print("A", xsh)
            oxsp = 1 + szx - xst
            xsp = szx
        else
    #         print("B", xsh)
            oxsp = 1 + xsp - xst + xsh
            # xsp = xsp
        end
        return xsp, oxsp
    end

    """
    Get cartesian indices for data from block cartesian indices.
    No out of bounds check is performed.
    """
    function data_sub_from_block_sub(block_size::ArrayOrTuple, margin_size::Int, bsub)
        xst = (block_size[1] * (bsub[1] - 1)) + 1 - margin_size
        xsp = (block_size[1] * (bsub[1] + 0)) + 0 + margin_size
        yst = (block_size[2] * (bsub[2] - 1)) + 1 - margin_size
        ysp = (block_size[2] * (bsub[2] + 0)) + 0 + margin_size
        zst = (block_size[3] * (bsub[3] - 1)) + 1 - margin_size
        zsp = (block_size[3] * (bsub[3] + 0)) + 0 + margin_size
        return xst, xsp, yst, ysp, zst, zsp
    end

    """
    Get face IDs based on position of cube in grid. Faces IDs along each axis
    are returned.

    f1, f2, f3 = get_face_ids_from_cube_in_grid([1,2,3], 1, false)
    """
    function get_face_ids_from_cube_in_grid(
        grid_size::ArrayOrTuple,
        cube_carthesian_position::ArrayOrTuple,
        trailing_face::Bool)
        if trailing_face
            trf = 1
        else
            trf = 0
        end
        sz1,sz2,sz3 = grid_size
        i, j, k = cube_carthesian_position
        f10 = (sz2 * sz3) * (i - 1 + trf)  + (j - 1) * sz3 + k
        nax1 = (1 + sz1) * sz2 * sz3
        f20 = nax1 +
            (sz2 + 1) * sz3 * (i - 1)  + (j - 1 + trf) * sz3 + k

        nax2 = sz1 * (1 + sz2) * sz3
        # @debug ("2st axis faces: ", nax2, ", ")
        nax3_layer = (sz3 + 1) * sz2 * (i - 1)
        nax3_row = (j - 1) * (sz3 + 1)
        # @debug ("3st axis faces in one layer and in one row: ",
        #     nax3_layer, " ", nax3_row,  "\n")
        @debug "number of axis1 faces: $nax1 axis2: $nax2 "
        f30 = nax1 +  nax2 +
            nax3_layer + nax3_row  + k + trf

        return f10, f20, f30

    end

    function cartesian_withloops(x,y)
        leny=length(y)
        lenx=length(x)
        m=leny*lenx
        OUT = zeros(eltype(x), m,2)
        c=1
        for i = 1:lenx
            for j = 1:leny
                OUT[c,1] = x[i]
                OUT[c,2] = y[j]
                c+=1
            end
        end
        return OUT
    end

    """
    Get list of surface faces in inner subgrid.
    """
    function cube_in_block_surface(
        block_size::ArrayOrTuple,
        cube_start::ArrayOrTuple,
        cube_stop::ArrayOrTuple
        )
        dimension = length(block_size)
        dim = 1
        inner_block_size = cube_stop - cube_start + ones(size(cube_stop))
        number_of_facelets_per_dim = zeros(Int, dimension)
        for i=1:dimension
    #         print("inner_block_size ", inner_block_size, "\n")
            ones_size = copy(inner_block_size)
            facelet_size_on_this_dim = [ones_size[j] for j=1:length(ones_size) if j != i]
    #         print("facelet_size", facelet_size_on_this_dim, "\n")
            number_of_facelets_per_dim[i] = prod(facelet_size_on_this_dim)

    #         copy
        end
    #     print("number_of_facelets_per_dim ", number_of_facelets_per_dim, "\n")
        total_number_of_facelets = sum(number_of_facelets_per_dim) * 2

        # output array
        facelet_inds = Array{Int64}(undef, total_number_of_facelets)

    #     Array
    #     array = Array{Int64}(undef, 5)
        ranges = [collect(cube_start[i]:cube_stop[i]) for i=1:dimension]
    #     print("ranges ", ranges, "\n")
    #     ranges = Array{Any}(undef, dimension)
        cart_index = zeros(Int64, dimension)
        linear_facelet_index = 1
        for i=1:dimension
            rest_dims = [j for j=1:dimension if j != i]
    #         print("rest dims ", rest_dims, "\n")
            r1 = collect(cube_start[rest_dims[1]]:cube_stop[rest_dims[1]])
            r2 = collect(cube_start[rest_dims[2]]:cube_stop[rest_dims[2]])
    #         print("r1 ", r1, ",", typeof(r1),"\n")
    #         print("r2 ", r2, ",", typeof(r2),"\n")
            cartrange_i_dim = cartesian_withloops( r1, r2, )
    #         println("cartrange_i_dim ", cartrange_i_dim)
            for k=1:size(cartrange_i_dim)[1]
                cart_index_rest = cartrange_i_dim[k, :]
                cart_index[i] = cube_start[i]
                for j=1:(dimension - 1)
                    cart_index[rest_dims[j]] = cart_index_rest[j]
                end
                facelet_inds[linear_facelet_index] = get_face_ids_from_cube_in_grid(
                    block_size, cart_index, false)[i]
                linear_facelet_index += 1
                cart_index[i] = cube_stop[i]
                facelet_inds[linear_facelet_index] = get_face_ids_from_cube_in_grid(
                    block_size, cart_index, true)[i]
    #             println("facelet_index: ", linear_facelet_index, " cart_index ", cart_index, " index rest ", cart_index_rest)
                linear_facelet_index += 1

            end

        end
        return facelet_inds
    end
    #     rest_dims = [2, 3]
    #     # shape of flat area (with dimension D-1). One dimension is 1 the other are
    #     # the same
    #     ones_size = copy(block_size)
    #     ones_size[dim] = 1
    #     # ones_size = ones(Int64, dimension)
    #     # ones_size[dim] = block_size[dim]
    #     ones(typeof(cube_start), ones_size...)
    #     # cat(dim, )
    #     # cartesian_withloops(collect(cubes_start[1]), collect(cube_stop[2]))
    #     # for
    #     return 1
    # end
    #
    # function a(i)
    #     return 1
    function grid_voxel_cart_to_id(grid_size::ArrayOrTuple, voxel_cart::ArrayOrTuple)
        sz1,sz2,sz3 = grid_size
        i, j, k = voxel_cart
        trf = 0
        ind = (sz2 * sz3) * (i - 1 + trf)  + (j - 1) * sz3 + k
        return ind
    end

    function grid_voxel_cart_to_node_id(grid_size::Array{Int, 1}, voxel_cart::ArrayOrTuple)
        sz1,sz2,sz3 = grid_size + [1,1,1]
        i, j, k = voxel_cart
        trf = 0
        ind = (sz2 * sz3) * (i - 1 + trf)  + (j - 1) * sz3 + k
        return ind

    end

    function grid_voxel_id_to_cart(grid_size::Array{Int, 1}, ind::Int64)
    #     ind = (sz2 * sz3) * (i - 1)  + (j - 1) * sz3 + k
        sz1,sz2,sz3 = grid_size
        layer = sz2*sz3
        ir = div(ind - 1, layer) + 1
        rest1 = mod(ind - 1 , layer) + 1
    #     println("rest1 ", rest1)
        row = sz3
        jr = div(rest1 -1, row) + 1
        rest2 = mod(rest1 - 1, row) + 1
        kr = rest2
    #     println("rests: ", rest1, " ", rest2)
        return [ir, jr, kr]
    end

    function grid_x_face_to_carthesian(grid_size::ArrayOrTuple, ind::Integer)
        return grid_voxel_id_to_cart(grid_size, ind)
    end

    function grid_y_face_to_carthesian(grid_size::ArrayOrTuple, ind::Integer)
    #     f20 = nax1 +
    #         (sz2 + 1) * sz3 * (i - 1 + trf)  + (j - 1 + trf) * sz3 + k
        sz1,sz2,sz3 = grid_size
        nax1 = (1 + sz1) * sz2 * sz3
        ind = ind - nax1

        layer = (sz2 + 1) * sz3
        ir = div(ind - 1, layer) + 1
        rest1 = mod(ind - 1, layer) + 1

        jr = div(rest1 - 1, sz3) + 1
        rest2 = mod(rest1 - 1, sz3) + 1
        kr = rest2
        return [ir, jr, kr]
    end

    function grid_z_face_to_carthesian(grid_size::ArrayOrTuple, ind::Int)
    #     nax3_layer = (sz3 + 1) * sz2 * (i - 1)
    #     nax3_row = (j - 1) * (sz3 + 1)
    #     f30 = nax1 +  nax2 +
    #         nax3_layer + nax3_row  + k + trf
        sz1,sz2,sz3 = grid_size
        nax1 = (1 + sz1) * sz2 * sz3
        nax2 = sz1 * (1 + sz2) * sz3
        ind = ind - nax1 - nax2

        layer = (sz3 + 1) * sz2
        ir = div(ind - 1, layer) + 1
        rest1 = mod(ind - 1, layer) + 1

        row = (sz3 + 1)
        jr = div(rest1 - 1, row) + 1
        rest2 = mod(rest1 - 1, row) + 1
        kr = rest2
        return [ir, jr, kr]
    end
    # end

    function grid_face_id_to_cartesian(grid_size::ArrayOrTuple, fid::Integer)
        sz1,sz2,sz3 = grid_size
        nax1 = (1 + sz1) * sz2 * sz3
        nax2 = sz1 * (1 + sz2) * sz3

        if fid <= nax1
            # x-face it is the same index as voxel index
            axis = 1
            voxel_cart = grid_x_face_to_carthesian(grid_size, fid)

        elseif fid <= (nax1 + nax2)
            # y-face
            axis = 2
            voxel_cart = grid_y_face_to_carthesian(grid_size, fid)
        else
            # z-face
            axis = 3
            voxel_cart = grid_z_face_to_carthesian(grid_size, fid)
        end
        return voxel_cart, axis
    end

    function grid_to_linear(data3d::AbstractArray, threshold=0)
        """
        Get grid linearized version of segmentation
        """
        segClin = spzeros(Int8, prod(size(data3d)), 1)
        sz = size(data3d)
        for k in 1:sz[3]
            for j in 1:sz[2]
                for i in 1:sz[1]
                    if data3d[i,j,k] > threshold
                        ind = grid_voxel_cart_to_id(size(data3d), [i,j,k])
                        segClin[ind] = 1
                    end
                end
            end
        end
        return segClin
    end

    """
    Get face ID in original grid from its ID in subgrid.
    """
    function sub_grid_face_id_to_orig_grid_face_id(
        data_size::ArrayOrTuple,
        block_size::ArrayOrTuple,
        offset::ArrayOrTuple,
        fid::Integer
        )
        face_cart, axis = grid_face_id_to_cartesian(block_size, fid)

        voxel_cart = face_cart + offset
        big_fids = get_face_ids_from_cube_in_grid(data_size, voxel_cart, false)
        big_fid = big_fids[axis]
        # print("face cart: $face_cart, axis $axis, new face ID: $big_fid, data size: $data_size")
        return big_fid, voxel_cart
    end

    """
    Calculate number of faces in grid.
    > grid_number_of_faces([5,5,5])

    """
    function grid_number_of_faces(grid_size::ArrayOrTuple)
        pr = prod(grid_size)
        num = pr + grid_size[1] * grid_size[2] +
            pr + grid_size[2] * grid_size[3] +
            pr + grid_size[1] * grid_size[3]
        return num

    end

    """
    Compute number of vertexes in grid.
    """
    function grid_number_of_nodes(grid_size::ArrayOrTuple)
        return (grid_size[1] + 1) * (grid_size[2] + 1) * (grid_size[3] + 1)
    end

    """
    Calculate node IDs from face ID in grid.
    Return array of node IDs and array of nodes carts.

    julia> nodes_ids, nodes_carts = grid_face_id_to_node_ids([2,3,4], 20)
    ([29, 34, 35, 30], Array{Int64,1}[[2, 2, 4], [2, 3, 4], [2, 2, 5], [2, 3, 5]])
    """
    function grid_face_id_to_node_ids(grid_size::ArrayOrTuple, face_id::Integer)
        voxel_cart, axis = grid_face_id_to_cartesian(grid_size, face_id)
        # println("voxel_cart, axis: ", voxel_cart, axis)
        axises = [1,2,3]
        deleteat!(axises, axis)
        voxel_cart1 = copy(voxel_cart)
        voxel_cart2 = copy(voxel_cart)
        voxel_cart3 = copy(voxel_cart)
        voxel_cart4 = copy(voxel_cart)

        voxel_cart2[axises[1]] += 1
        voxel_cart3[axises[1]] += 1
        voxel_cart3[axises[2]] += 1
        voxel_cart4[axises[2]] += 1

        ndi1 = grid_voxel_cart_to_node_id(grid_size, voxel_cart1)
        ndi2 = grid_voxel_cart_to_node_id(grid_size, voxel_cart2)
        ndi3 = grid_voxel_cart_to_node_id(grid_size, voxel_cart3)
        ndi4 = grid_voxel_cart_to_node_id(grid_size, voxel_cart4)
        carts = [voxel_cart1, voxel_cart2, voxel_cart3, voxel_cart4]
        return [ndi1, ndi2, ndi3, ndi4], carts
        # TODO finish

    end

    function grid_Fchar_to_V_FVreduced(Fchar::SparseArrays.SparseVector, data_size::Tuple)
        data_size = size_as_array(data_size)
        return grid_Fchar_to_V_FVreduced(Fchar, data_size)
    end

    """
    Calculate V and FV based on linear characteristic matrix of F and size of data.
    """
    function grid_Fchar_to_V_FVreduced(Fchar::SparseArrays.SparseVector, data_size::Array)
        all_info = [
            grid_face_id_to_node_ids(data_size, i)
            for i=1:length(Fchar) if Fchar[i] == 1
        ]

        filtered_bigFV2 = [all_info[i][1] for i=1:length(all_info)]
        n_nodes = grid_number_of_nodes(data_size)
        # Vcomputed = spzeros(Float64, 3, n_nodes)
        Vcomputed = zeros(Float64, 3, n_nodes)
        for i=1:length(all_info)
            # all_info[i][2]
            Vcomputed[:, all_info[i][1][1]] = all_info[i][2][1]
            Vcomputed[:, all_info[i][1][2]] = all_info[i][2][2]
            Vcomputed[:, all_info[i][1][3]] = all_info[i][2][3]
            Vcomputed[:, all_info[i][1][4]] = all_info[i][2][4]
        end
        return Vcomputed, filtered_bigFV2
    end

    function count_F_from_Fchar(Fchar)
        countF = 0
        for i=1:length(Fchar)
            if Fchar[i] == 1
                countF
                countF = countF + 1
            end
        end
        return countF
    end

    """
    Calculate V and FV based on linear characteristic matrix of F and size of data.
    """
    function grid_Fchar_to_Vreduced_FVreduced(Fchar::SparseArrays.SparseVector, data_size::Array;
        voxelsize::Array=[1,1,1]
        )
        # data_size = size_as_array(size(segmentation))
        countF = nnz(Fchar)

        # countF = count_F_from_Fchar(Fchar)

        # Construct FV with indexes to full V
        FV3 = Array{Array{Int64,1},1}(undef, countF)
        node_carts_dict = Dict()
        fv_i = 0
        rows, vals = findnz(Fchar)
        # println("=== findnz ")
        # display(Fchar)
        for i in rows
            # print("$i ")
            face_ids, nodes_carts = grid_face_id_to_node_ids(data_size, i)
            node_carts_dict[face_ids[1]] = nodes_carts[1]
            node_carts_dict[face_ids[2]] = nodes_carts[2]
            node_carts_dict[face_ids[3]] = nodes_carts[3]
            node_carts_dict[face_ids[4]] = nodes_carts[4]
            fv_i += 1
            FV3[fv_i] = face_ids
        end

        # for i=1:length(Fchar)
        #     if Fchar[i] == 1
        #         face_ids, nodes_carts = grid_face_id_to_node_ids(data_size, i)
        #         node_carts_dict[face_ids[1]] = nodes_carts[1]
        #         node_carts_dict[face_ids[2]] = nodes_carts[2]
        #         node_carts_dict[face_ids[3]] = nodes_carts[3]
        #         node_carts_dict[face_ids[4]] = nodes_carts[4]
        #         fv_i += 1
        #         FV3[fv_i] = face_ids
        #     end
        # end
        # create V reduced
        V3arr = collect(values(node_carts_dict))
        V3 = Array{Float64,2}(undef, 3, length(V3arr))
        for i=1:length(V3arr)
            V3[:, i] = V3arr[i] .* voxelsize
        end

        # Relabel FV with full V to FV with reduced V
        ks = keys(node_carts_dict)
        node_ids = Dict(zip(ks, collect(1:length(ks))))

        # convert
        for i=1:length(FV3)

            for j=1:length(FV3[i])
                FV3[i][j] = node_ids[FV3[i][j]]
            end
        end
        return V3, FV3
    end

    function grid_Fchar_to_Vreduced_FVreduced_old(Fchar::SparseArrays.SparseVector, data_size::Array)
        # data_size = size_as_array(size(segmentation))

        countF = count_F_from_Fchar(Fchar)

        # Construct FV with indexes to full V
        FV3 = Array{Array{Int64,1},1}(undef, countF)
        node_carts_dict = Dict()
        fv_i = 0
        for i=1:length(Fchar)
            if Fchar[i] == 1
                face_ids, nodes_carts = grid_face_id_to_node_ids(data_size, i)
                node_carts_dict[face_ids[1]] = nodes_carts[1]
                node_carts_dict[face_ids[2]] = nodes_carts[2]
                node_carts_dict[face_ids[3]] = nodes_carts[3]
                node_carts_dict[face_ids[4]] = nodes_carts[4]
                fv_i += 1
                FV3[fv_i] = face_ids
            end
        end
        # create V reduced
        V3arr = collect(values(node_carts_dict))
        V3 = Array{Float64,2}(undef, 3, length(V3arr))
        for i=1:length(V3arr)
            V3[:, i] = V3arr[i]
        end

        # Relabel FV with full V to FV with reduced V
        ks = keys(node_carts_dict)
        node_ids = Dict(zip(ks, collect(1:length(ks))))

        # convert
        for i=1:length(FV3)

            for j=1:length(FV3[i])
                FV3[i][j] = node_ids[FV3[i][j]]
            end
        end
        return V3, FV3
    end

    """
    Calculate V and FV based on linear characteristic matrix of F and size of data.
    It is calculated full FV and than it is reduced
    """
    function grid_Fchar_to_V_FVfulltoreduced(Fchar::SparseArrays.SparseVector, data_size::Array)
        bigV, topology = Lar.cuboidGrid(data_size, true)
        (bigVV, bigEV, bigFV, bigCV) = topology

        __fix_cuboidGrid_FV!((bigV, topology))

        # Get FV and filter double faces on the border
        filtered_bigFV = [
            bigFV[i] for i=1:length(Fchar) if Fchar[i] == 1
        ]
        return bigV, filtered_bigFV

    end

# end


text:
# First Half

filter_lf(s) = filter(c -> c ≠ '\n', collect(s))

sumcounts1(strs) = map(length ∘ unique ∘ sort ∘ filter_lf, strs) |> sum
# For each group, count the number of questions ANYONE answered "yes"
# and then sum those counts

puzzleinput = open("inputs/day06.txt") do f
    split(read(f, String), "\n\n")
end

sumcounts1(puzzleinput)

# Second Half

spliton_lf(s) = split(s, "\n") |> collect
intersect_qs(cs) = foldl(intersect, cs, init=collect('a':'z'))

sumcounts2(strs) = map(length ∘ intersect_qs ∘ spliton_lf, strs) |> sum
# For each group, count the number of questions EVERYONE answered "yes"
# and then sum those counts

puzzleinput = open("inputs/day06.txt") do f
    split(read(f, String), "\n\n")
end

sumcounts2(puzzleinput)


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

mutable struct IoK8sApiCoreV1PersistentVolume <: SwaggerModel
    apiVersion::Any # spec type: Union{ Nothing, String } # spec name: apiVersion
    kind::Any # spec type: Union{ Nothing, String } # spec name: kind
    metadata::Any # spec type: Union{ Nothing, IoK8sApimachineryPkgApisMetaV1ObjectMeta } # spec name: metadata
    spec::Any # spec type: Union{ Nothing, IoK8sApiCoreV1PersistentVolumeSpec } # spec name: spec
    status::Any # spec type: Union{ Nothing, IoK8sApiCoreV1PersistentVolumeStatus } # spec name: status

    function IoK8sApiCoreV1PersistentVolume(;apiVersion=nothing, kind=nothing, metadata=nothing, spec=nothing, status=nothing)
        o = new()
        validate_property(IoK8sApiCoreV1PersistentVolume, Symbol("apiVersion"), apiVersion)
        setfield!(o, Symbol("apiVersion"), apiVersion)
        validate_property(IoK8sApiCoreV1PersistentVolume, Symbol("kind"), kind)
        setfield!(o, Symbol("kind"), kind)
        validate_property(IoK8sApiCoreV1PersistentVolume, Symbol("metadata"), metadata)
        setfield!(o, Symbol("metadata"), metadata)
        validate_property(IoK8sApiCoreV1PersistentVolume, Symbol("spec"), spec)
        setfield!(o, Symbol("spec"), spec)
        validate_property(IoK8sApiCoreV1PersistentVolume, Symbol("status"), status)
        setfield!(o, Symbol("status"), status)
        o
    end
end # type IoK8sApiCoreV1PersistentVolume

const _property_map_IoK8sApiCoreV1PersistentVolume = Dict{Symbol,Symbol}(Symbol("apiVersion")=>Symbol("apiVersion"), Symbol("kind")=>Symbol("kind"), Symbol("metadata")=>Symbol("metadata"), Symbol("spec")=>Symbol("spec"), Symbol("status")=>Symbol("status"))
const _property_types_IoK8sApiCoreV1PersistentVolume = Dict{Symbol,String}(Symbol("apiVersion")=>"String", Symbol("kind")=>"String", Symbol("metadata")=>"IoK8sApimachineryPkgApisMetaV1ObjectMeta", Symbol("spec")=>"IoK8sApiCoreV1PersistentVolumeSpec", Symbol("status")=>"IoK8sApiCoreV1PersistentVolumeStatus")
Base.propertynames(::Type{ IoK8sApiCoreV1PersistentVolume }) = collect(keys(_property_map_IoK8sApiCoreV1PersistentVolume))
Swagger.property_type(::Type{ IoK8sApiCoreV1PersistentVolume }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_IoK8sApiCoreV1PersistentVolume[name]))}
Swagger.field_name(::Type{ IoK8sApiCoreV1PersistentVolume }, property_name::Symbol) =  _property_map_IoK8sApiCoreV1PersistentVolume[property_name]

function check_required(o::IoK8sApiCoreV1PersistentVolume)
    true
end

function validate_property(::Type{ IoK8sApiCoreV1PersistentVolume }, name::Symbol, val)
end


text:
int_rules_1_1_4_3 = @theory begin

    #= ::Subsection::Closed:: =#
    #= 1.1.4.3*(e*x)^m*(a*x^j+b*x^k)^p*(c+d*x^n)^q =#
    @apply_utils Antiderivative((~x) ^ ~(m') * (~(a') * (~x) ^ ~j + ~(b') * (~x) ^ ~(k')) ^ ~p * (~c + ~(d') * (~x) ^ ~n) ^ ~(q'), ~x) => (1 / ~n) * Subst(Antiderivative((~x) ^ (Simplify((~m + 1) / ~n) - 1) * (~a * (~x) ^ Simplify(~j / ~n) + ~b * (~x) ^ Simplify(~k / ~n)) ^ ~p * (~c + ~d * ~x) ^ ~q, ~x), ~x, (~x) ^ ~n) <-- FreeQ([~a, ~b, ~c, ~d, ~j, ~k, ~m, ~n, ~p, ~q], ~x) && (Not(IntegerQ(~p)) && (NeQ(~k, ~j) && (IntegerQ(Simplify(~j / ~n)) && (IntegerQ(Simplify(~k / ~n)) && (IntegerQ(Simplify((~m + 1) / ~n)) && NeQ((~n) ^ 2, 1))))))
    @apply_utils Antiderivative((~e * ~x) ^ ~(m') * (~(a') * (~x) ^ ~j + ~(b') * (~x) ^ ~(k')) ^ ~p * (~c + ~(d') * (~x) ^ ~(n')) ^ ~(q'), ~x) => (((~e) ^ IntPart(~m) * (~e * ~x) ^ FracPart(~m)) / (~x) ^ FracPart(~m)) * Antiderivative((~x) ^ ~m * (~a * (~x) ^ ~j + ~b * (~x) ^ ~k) ^ ~p * (~c + ~d * (~x) ^ ~n) ^ ~q, ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~j, ~k, ~m, ~n, ~p, ~q], ~x) && (Not(IntegerQ(~p)) && (NeQ(~k, ~j) && (IntegerQ(Simplify(~j / ~n)) && (IntegerQ(Simplify(~k / ~n)) && (IntegerQ(Simplify((~m + 1) / ~n)) && NeQ((~n) ^ 2, 1))))))
    @apply_utils Antiderivative((~(e') * ~x) ^ ~(m') * (~(a') * (~x) ^ ~(j') + ~(b') * (~x) ^ ~(jn')) ^ ~p * (~c + ~(d') * (~x) ^ ~(n')), ~x) => (~c * (~e) ^ (~j - 1) * (~e * ~x) ^ ((~m - ~j) + 1) * (~a * (~x) ^ ~j + ~b * (~x) ^ (~j + ~n)) ^ (~p + 1)) / (~a * (~m + ~j * ~p + 1)) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~j, ~m, ~n, ~p], ~x) && (EqQ(~jn, ~j + ~n) && (Not(IntegerQ(~p)) && (NeQ(~b * ~c - ~a * ~d, 0) && (EqQ(~a * ~d * (~m + ~j * ~p + 1) - ~b * ~c * (~m + ~n + ~p * (~j + ~n) + 1), 0) && ((GtQ(~e, 0) || IntegersQ(~j)) && NeQ(~m + ~j * ~p + 1, 0))))))
    @apply_utils Antiderivative((~(e') * ~x) ^ ~(m') * (~(a') * (~x) ^ ~(j') + ~(b') * (~x) ^ ~(jn')) ^ ~p * (~c + ~(d') * (~x) ^ ~(n')), ~x) => (-((~e) ^ (~j - 1)) * (~b * ~c - ~a * ~d) * (~e * ~x) ^ ((~m - ~j) + 1) * (~a * (~x) ^ ~j + ~b * (~x) ^ (~j + ~n)) ^ (~p + 1)) / (~a * ~b * ~n * (~p + 1)) - (((~e) ^ ~j * (~a * ~d * (~m + ~j * ~p + 1) - ~b * ~c * (~m + ~n + ~p * (~j + ~n) + 1))) / (~a * ~b * ~n * (~p + 1))) * Antiderivative((~e * ~x) ^ (~m - ~j) * (~a * (~x) ^ ~j + ~b * (~x) ^ (~j + ~n)) ^ (~p + 1), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~j, ~m, ~n], ~x) && (EqQ(~jn, ~j + ~n) && (Not(IntegerQ(~p)) && (NeQ(~b * ~c - ~a * ~d, 0) && (LtQ(~p, -1) && (GtQ(~j, 0) && (LeQ(~j, ~m) && (GtQ(~e, 0) || IntegerQ(~j))))))))
    @apply_utils Antiderivative((~(e') * ~x) ^ ~(m') * (~(a') * (~x) ^ ~(j') + ~(b') * (~x) ^ ~(jn')) ^ ~p * (~c + ~(d') * (~x) ^ ~(n')), ~x) => (~c * (~e) ^ (~j - 1) * (~e * ~x) ^ ((~m - ~j) + 1) * (~a * (~x) ^ ~j + ~b * (~x) ^ (~j + ~n)) ^ (~p + 1)) / (~a * (~m + ~j * ~p + 1)) + ((~a * ~d * (~m + ~j * ~p + 1) - ~b * ~c * (~m + ~n + ~p * (~j + ~n) + 1)) / (~a * (~e) ^ ~n * (~m + ~j * ~p + 1))) * Antiderivative((~e * ~x) ^ (~m + ~n) * (~a * (~x) ^ ~j + ~b * (~x) ^ (~j + ~n)) ^ ~p, ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~j, ~p], ~x) && (EqQ(~jn, ~j + ~n) && (Not(IntegerQ(~p)) && (NeQ(~b * ~c - ~a * ~d, 0) && (GtQ(~n, 0) && ((LtQ(~m + ~j * ~p, -1) || IntegersQ(~m - 1 / 2, ~p - 1 / 2) && (LtQ(~p, 0) && LtQ(~m, -(~n) * ~p - 1))) && ((GtQ(~e, 0) || IntegersQ(~j, ~n)) && (NeQ(~m + ~j * ~p + 1, 0) && NeQ((~m - ~n) + ~j * ~p + 1, 0))))))))
    @apply_utils Antiderivative((~(e') * ~x) ^ ~(m') * (~(a') * (~x) ^ ~(j') + ~(b') * (~x) ^ ~(jn')) ^ ~p * (~c + ~(d') * (~x) ^ ~(n')), ~x) => (~d * (~e) ^ (~j - 1) * (~e * ~x) ^ ((~m - ~j) + 1) * (~a * (~x) ^ ~j + ~b * (~x) ^ (~j + ~n)) ^ (~p + 1)) / (~b * (~m + ~n + ~p * (~j + ~n) + 1)) - ((~a * ~d * (~m + ~j * ~p + 1) - ~b * ~c * (~m + ~n + ~p * (~j + ~n) + 1)) / (~b * (~m + ~n + ~p * (~j + ~n) + 1))) * Antiderivative((~e * ~x) ^ ~m * (~a * (~x) ^ ~j + ~b * (~x) ^ (~j + ~n)) ^ ~p, ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~j, ~m, ~n, ~p], ~x) && (EqQ(~jn, ~j + ~n) && (Not(IntegerQ(~p)) && (NeQ(~b * ~c - ~a * ~d, 0) && (NeQ(~m + ~n + ~p * (~j + ~n) + 1, 0) && (GtQ(~e, 0) || IntegerQ(~j))))))
    @apply_utils Antiderivative((~x) ^ ~(m') * (~(a') * (~x) ^ ~j + ~(b') * (~x) ^ ~(k')) ^ ~p * (~c + ~(d') * (~x) ^ ~(n')) ^ ~(q'), ~x) => (1 / (~m + 1)) * Subst(Antiderivative((~a * (~x) ^ Simplify(~j / (~m + 1)) + ~b * (~x) ^ Simplify(~k / (~m + 1))) ^ ~p * (~c + ~d * (~x) ^ Simplify(~n / (~m + 1))) ^ ~q, ~x), ~x, (~x) ^ (~m + 1)) <-- FreeQ([~a, ~b, ~c, ~d, ~j, ~k, ~m, ~n, ~p, ~q], ~x) && (Not(IntegerQ(~p)) && (NeQ(~k, ~j) && (IntegerQ(Simplify(~j / ~n)) && (IntegerQ(Simplify(~k / ~n)) && (NeQ(~m, -1) && (IntegerQ(Simplify(~n / (~m + 1))) && Not(IntegerQ(~n))))))))
    @apply_utils Antiderivative((~e * ~x) ^ ~(m') * (~(a') * (~x) ^ ~j + ~(b') * (~x) ^ ~(k')) ^ ~p * (~c + ~(d') * (~x) ^ ~(n')) ^ ~(q'), ~x) => (((~e) ^ IntPart(~m) * (~e * ~x) ^ FracPart(~m)) / (~x) ^ FracPart(~m)) * Antiderivative((~x) ^ ~m * (~a * (~x) ^ ~j + ~b * (~x) ^ ~k) ^ ~p * (~c + ~d * (~x) ^ ~n) ^ ~q, ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~j, ~k, ~m, ~n, ~p, ~q], ~x) && (Not(IntegerQ(~p)) && (NeQ(~k, ~j) && (IntegerQ(Simplify(~j / ~n)) && (IntegerQ(Simplify(~k / ~n)) && (NeQ(~m, -1) && (IntegerQ(Simplify(~n / (~m + 1))) && Not(IntegerQ(~n))))))))
    @apply_utils Antiderivative((~(e') * ~x) ^ ~(m') * (~(a') * (~x) ^ ~(j') + ~(b') * (~x) ^ ~(jn')) ^ ~p * (~c + ~(d') * (~x) ^ ~(n')) ^ ~(q'), ~x) => (((~e) ^ IntPart(~m) * (~e * ~x) ^ FracPart(~m) * (~a * (~x) ^ ~j + ~b * (~x) ^ (~j + ~n)) ^ FracPart(~p)) / ((~x) ^ (FracPart(~m) + ~j * FracPart(~p)) * (~a + ~b * (~x) ^ ~n) ^ FracPart(~p))) * Antiderivative((~x) ^ (~m + ~j * ~p) * (~a + ~b * (~x) ^ ~n) ^ ~p * (~c + ~d * (~x) ^ ~n) ^ ~q, ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~j, ~m, ~n, ~p, ~q], ~x) && (EqQ(~jn, ~j + ~n) && (Not(IntegerQ(~p)) && (NeQ(~b * ~c - ~a * ~d, 0) && Not(EqQ(~n, 1) && EqQ(~j, 1)))))
end

max_stars_repo_path: src/rules/1 Algebraic functions/1.1 Binomial products/1.1.4 Improper/1.1.4.3 (e x)^m (a x^j+b x^k)^p (c+d x^n)^q.jl


text:
module Examples

using ..Common

function slab_pml(; D::Real=0, ω::Real = 40, γ::Real=2, n::Real=3, dx::Real=2π/ω/n/20)
    cavity = NondispersiveDomain(Interval(-.5,.5),n)
    pump = DispersiveDomain(Interval(-.5,.5),TwoLevelSystem(ω,D,γ))
    bnd = Boundary(Interval(-1,1),PML,DirichletBC)
    lat = Cartesian(dx)
    latdom = LatticeDomain(bnd,lat)
    return smooth!(Simulation(ω,cavity,pump,latdom))
end

function slab_cpml(; D::Real=0, ω::Real = 40, γ::Real=2, n::Real=3, dx::Real=2π/ω/n/20)
    cavity = NondispersiveDomain(Interval(-.5,.5),n)
    pump = DispersiveDomain(Interval(-.5,.5),TwoLevelSystem(ω,D,γ))
    bnd = Boundary(Interval(-1,1),cPML,DirichletBC)
    lat = Cartesian(dx)
    latdom = LatticeDomain(bnd,lat)
    return smooth!(Simulation(ω,cavity,pump,latdom))
end

function slab_out(; D::Real=0, ω::Real = 40, γ::Real=2, n::Real=3, dx::Real=2π/ω/n/20)
    cavity = NondispersiveDomain(Interval(-.5,.5),n)
    pump = DispersiveDomain(Interval(-.5,.5),TwoLevelSystem(ω,D,γ))
    bnd = Boundary(Interval(-1,1),MatchedBC{1}(out=[1]),MatchedBC{2}(out=[1]))
    lat = Cartesian(dx)
    latdom = LatticeDomain(bnd,lat)
    return smooth!(Simulation(ω,cavity,pump,latdom))
end

function slab_in(; D::Real=0, ω::Real = 40, γ::Real=2, n::Real=3, dx::Real=2π/ω/n/20)
    cavity = NondispersiveDomain(Interval(-.5,.5),n)
    pump = DispersiveDomain(Interval(-.5,.5),TwoLevelSystem(ω,D,γ))
    bnd = Boundary(Interval(-1,1),MatchedBC{1}(in=[1]),MatchedBC{2}(in=[1]))
    lat = Cartesian(dx)
    latdom = LatticeDomain(bnd,lat)
    return smooth!(Simulation(ω,cavity,pump,latdom))
end

end # module


text:
using Colors, ColorSchemes
using CairoMakie, FileIO
using CSV, DataFrames
using Statistics

imgpath = joinpath("Molenaar", "Figure 3", "Plots")
resultspath = joinpath("Molenaar", "Results")

x = [res["mu"] for res in results]
y = [res["used_ribo"] for res in results]

a = [] # transporter rate
b = [] # matabolic rate
c = [] # ribosomale rate
d = [] # lipid biosynthesis rate

for i = 1:100
    push!(a, y[i][1])
    push!(b, y[i][2])
    push!(c, y[i][3])
    push!(d, y[i][4])
end

a = Float64.(a)
b = Float64.(b)
c = Float64.(c)
d = Float64.(d)

metabolism = b
ribosomes = c
transporters = a
lipids = d

fracs = [metabolism, ribosomes, transporters, lipids]
fraclabels = ["Metabolism", "Ribosomes", "Transporters", "Lipidsynthesis"]
function plot_proteome(ax, x, fracs, fraclabels)
    d = cumsum(hcat(fracs...)'; dims = 1)
    d = [zeros(size(d, 2))'; d ./ d[end, :]']
    for i = 1:length(fraclabels)
        band!(ax, x, d[i, :], d[i+1, :], label = fraclabels[i])
    end
end

f = Figure();
ax = Axis(f[1, 1]);
plot_proteome(ax, 1:100, fracs, fraclabels)
ax.xlabel = "Simulations with rising substarte concentration"
ax.ylabel = "Proteome fraction"
f[1, 2] = Legend(f, ax, "Protein", unique = true, framevisible = false)
f

FileIO.save(joinpath(imgpath, "Figure_3Cband.pdf"), f)


text:
abstract type AbstractDiffusion <: AbstractPhysicsComponent end

struct ConstantViscosity{𝒯} <: AbstractTerm{𝒯} end

@inline function calc_diffusive_flux_argument!(grad, ::Nothing, _...)
    grad.∇ρ = 0
    grad.∇u = @SVector [0, 0, 0]
    grad.∇θ = 0

    return nothing
end

@inline function calc_diffusive_flux_argument!(grad, diff::ConstantViscosity, state, aux, physics)
    ρ = state.ρ
    ρu = state.ρu
    ρθ = state.ρθ

    u = ρu / ρ
    θ = ρθ / ρ

    grad.∇ρ = ρ
    grad.∇u = u
    grad.∇θ = θ

    return nothing
end

@inline function calc_diffusive_flux!(gradflux, ::Nothing, _...)
    gradflux.μ∇ρ = @SVector [0, 0, 0]
    gradflux.ν∇u = @SMatrix zeros(3, 3)
    gradflux.κ∇θ = @SVector [0, 0, 0]

    return nothing
end

@inline function calc_diffusive_flux!(gradflux, ::ConstantViscosity, grad, state, aux, physics)
    μ = physics.params.μ * I
    ν = physics.params.ν * I
    κ = physics.params.κ * I

    gradflux.μ∇ρ = -μ * grad.∇ρ
    gradflux.ν∇u = -ν * grad.∇u
    gradflux.κ∇θ = -κ * grad.∇θ

    return nothing
end

max_stars_repo_name: CliMA/CouplerMachine

max_issues_repo_name: CliMA/CouplerMachine


text:
function create_stages(template_uc, c_sys5_uc, c_sys5_ed)
    return Dict(
        "UC" => Stage(GenericOpProblem, template_uc, c_sys5_uc, GLPK_optimizer),
        "ED" => Stage(GenericOpProblem, template_ed, c_sys5_ed, GLPK_optimizer),
    )
end

function create_sequence()
    return SimulationSequence(
        step_resolution = Hour(24),
        order = Dict(1 => "UC", 2 => "ED"),
        feedforward_chronologies = Dict(("UC" => "ED") => Synchronize(periods = 24)),
        horizons = Dict("UC" => 24, "ED" => 12),
        intervals = Dict(
            "UC" => (Hour(24), Consecutive()),
            "ED" => (Hour(1), Consecutive()),
        ),
        feedforward = Dict(
            ("ED", :devices, :Generators) => SemiContinuousFF(
                binary_source_stage = PSI.ON,
                affected_variables = [PSI.ACTIVE_POWER],
            ),
        ),
        ini_cond_chronology = InterStageChronology(),
    )
end

function test_sequence_build(file_path::String)
    c_sys5_uc = build_system("c_sys5_uc")
    c_sys5_ed = build_system("c_sys5_ed")

    @testset "Test Simulation Simulation Sequence Validation" begin
        sequence = create_sequence()
        @test length(findall(x -> x == 2, sequence.execution_order)) == 24
        @test length(findall(x -> x == 1, sequence.execution_order)) == 1
    end

    @testset "Simulation with provided initial time" begin
        stages_definition = create_stages(template_basic_uc, c_sys5_uc, c_sys5_ed)
        sequence = create_sequence()
        second_day = DayAhead[24] + Hour(1)
        sim = Simulation(
            name = "test",
            steps = 1,
            stages = stages_definition,
            stages_sequence = sequence,
            simulation_folder = file_path,
            initial_time = second_day,
        )
        build!(sim)

        for stage in values(sim.stages)
            @test PSI.get_initial_time(stage) == second_day
        end
    end

    @testset "Simulation Build Tests" begin
        stages_definition = create_stages(template_basic_uc, c_sys5_uc, c_sys5_ed)
        sequence = create_sequence()
        sim = Simulation(
            name = "test",
            steps = 1,
            stages = stages_definition,
            stages_sequence = sequence,
            simulation_folder = file_path,
        )
        build!(sim)

        @test isempty(values(sim.internal.simulation_cache))
        for field in fieldnames(SimulationSequence)
            if fieldtype(SimulationSequence, field) == Union{Dates.DateTime, Nothing}
                @test !isnothing(getfield(sim.sequence, field))
            end
        end
        @test isa(sim.sequence, SimulationSequence)
    end

    ####################### Negative Tests ########################################
    @testset "Test when a simulation has incorrect arguments" begin
        @test_throws UndefKeywordError sim =
            Simulation(name = "test", steps = 1, simulation_folder = file_path)
    end

    @testset "Test if a wrong initial time is provided" begin
        stages_definition = create_stages(template_basic_uc, c_sys5_uc, c_sys5_ed)
        sequence = create_sequence()
        sim = Simulation(
            name = "test",
            steps = 1,
            stages = stages_definition,
            stages_sequence = sequence,
            simulation_folder = file_path,
            initial_time = Dates.now(),
        )
        @test_throws IS.ConflictingInputsError build!(sim)
    end

    @testset "Test if file path is not writeable" begin
        stages_definition = create_stages(template_basic_uc, c_sys5_uc, c_sys5_ed)
        sequence = SimulationSequence(
            step_resolution = Hour(24),
            order = Dict(1 => "UC", 2 => "ED"),
            feedforward_chronologies = Dict(("UC" => "ED") => Synchronize(periods = 24)),
            horizons = Dict("UC" => 24, "ED" => 12),
            intervals = Dict(
                "UC" => (Hour(24), Consecutive()),
                "ED" => (Hour(1), Consecutive()),
            ),
            feedforward = Dict(
                ("ED", :devices, :Generators) => SemiContinuousFF(
                    binary_source_stage = PSI.ON,
                    affected_variables = [PSI.ACTIVE_POWER],
                ),
            ),
            ini_cond_chronology = InterStageChronology(),
        )
        sim = Simulation(
            name = "fake_path",
            steps = 1,
            stages = stages_definition,
            stages_sequence = sequence,
            simulation_folder = "fake_path",
        )
        @test_throws IS.ConflictingInputsError PSI._check_folder(sim)
    end

    @testset "chronology look ahead length is too long for horizon" begin
        stages_definition = create_stages(template_basic_uc, c_sys5_uc, c_sys5_ed)
        sequence = SimulationSequence(
            step_resolution = Hour(24),
            order = Dict(1 => "UC", 2 => "ED"),
            feedforward_chronologies = Dict(("UC" => "ED") => Synchronize(periods = 30)),
            horizons = Dict("UC" => 24, "ED" => 12),
            intervals = Dict(
                "UC" => (Hour(24), Consecutive()),
                "ED" => (Hour(1), Consecutive()),
            ),
            feedforward = Dict(
                ("ED", :devices, :Generators) => SemiContinuousFF(
                    binary_source_stage = PSI.ON,
                    affected_variables = [PSI.ACTIVE_POWER],
                ),
            ),
            ini_cond_chronology = InterStageChronology(),
        )
        sim = Simulation(
            name = "look_ahead",
            steps = 1,
            stages = stages_definition,
            stages_sequence = sequence,
            simulation_folder = file_path,
        )
        @test_throws IS.ConflictingInputsError PSI._check_feedforward_chronologies(sim)
    end

    @testset "too long of a horizon for forecast" begin
        stages_definition = create_stages(template_basic_uc, c_sys5_uc, c_sys5_ed)
        sequence = SimulationSequence(
            step_resolution = Hour(24),
            order = Dict(1 => "UC", 2 => "ED"),
            feedforward_chronologies = Dict(("UC" => "ED") => Synchronize(periods = 24)),
            horizons = Dict("UC" => 72, "ED" => 12),
            intervals = Dict(
                "UC" => (Hour(24), Consecutive()),
                "ED" => (Hour(1), Consecutive()),
            ),
            feedforward = Dict(
                ("ED", :devices, :Generators) => SemiContinuousFF(
                    binary_source_stage = PSI.ON,
                    affected_variables = [PSI.ACTIVE_POWER],
                ),
            ),
            ini_cond_chronology = InterStageChronology(),
        )
        sim = Simulation(
            name = "long_horizon",
            steps = 1,
            stages = stages_definition,
            stages_sequence = sequence,
            simulation_folder = file_path,
        )
        sim.internal = PSI.SimulationInternal(
            sim.steps,
            keys(sim.sequence.order),
            mktempdir(),
            PSI.get_name(sim),
        )
        @test_throws IS.ConflictingInputsError PSI._get_simulation_initial_times!(sim)
    end

    @testset "Test too many steps for forecast" begin
        stages_definition = create_stages(template_basic_uc, c_sys5_uc, c_sys5_ed)
        sequence = SimulationSequence(
            step_resolution = Hour(24),
            order = Dict(1 => "UC", 2 => "ED"),
            feedforward_chronologies = Dict(("UC" => "ED") => Synchronize(periods = 24)),
            horizons = Dict("UC" => 24, "ED" => 12),
            intervals = Dict(
                "UC" => (Hour(24), Consecutive()),
                "ED" => (Hour(1), Consecutive()),
            ),
            feedforward = Dict(
                ("ED", :devices, :Generators) => SemiContinuousFF(
                    binary_source_stage = PSI.ON,
                    affected_variables = [PSI.ACTIVE_POWER],
                ),
            ),
            ini_cond_chronology = InterStageChronology(),
        )
        sim = Simulation(
            name = "steps",
            steps = 5,
            stages = stages_definition,
            stages_sequence = sequence,
            simulation_folder = file_path,
        )
        @test_throws IS.ConflictingInputsError build!(sim)
    end

    @testset "Test Creation of Simulations with Cache" begin
        stages_definition = create_stages(template_standard_uc, c_sys5_uc, c_sys5_ed)

        # Cache is not defined all together
        sequence_no_cache = SimulationSequence(
            step_resolution = Hour(24),
            order = Dict(1 => "UC", 2 => "ED"),
            feedforward_chronologies = Dict(("UC" => "ED") => Synchronize(periods = 24)),
            horizons = Dict("UC" => 24, "ED" => 12),
            intervals = Dict(
                "UC" => (Hour(24), Consecutive()),
                "ED" => (Hour(1), Consecutive()),
            ),
            feedforward = Dict(
                ("ED", :devices, :Generators) => SemiContinuousFF(
                    binary_source_stage = PSI.ON,
                    affected_variables = [PSI.ACTIVE_POWER],
                ),
            ),
            ini_cond_chronology = InterStageChronology(),
        )
        sim = Simulation(
            name = "cache",
            steps = 1,
            stages = stages_definition,
            stages_sequence = sequence_no_cache,
            simulation_folder = file_path,
        )
        build!(sim)

        @test !isempty(sim.internal.simulation_cache)

        stages_definition = create_stages(template_standard_uc, c_sys5_uc, c_sys5_ed)
        sequence = SimulationSequence(
            step_resolution = Hour(24),
            order = Dict(1 => "UC", 2 => "ED"),
            feedforward_chronologies = Dict(("UC" => "ED") => Synchronize(periods = 24)),
            horizons = Dict("UC" => 24, "ED" => 12),
            intervals = Dict(
                "UC" => (Hour(24), Consecutive()),
                "ED" => (Hour(1), Consecutive()),
            ),
            feedforward = Dict(
                ("ED", :devices, :Generators) => SemiContinuousFF(
                    binary_source_stage = PSI.ON,
                    affected_variables = [PSI.ACTIVE_POWER],
                ),
            ),
            cache = Dict(("UC",) => TimeStatusChange(PSY.ThermalStandard, PSI.ON)),
            ini_cond_chronology = InterStageChronology(),
        )
        sim = Simulation(
            name = "caches",
            steps = 2,
            stages = stages_definition,
            stages_sequence = sequence,
            simulation_folder = file_path,
        )

        build!(sim)

        @test !isempty(sim.internal.simulation_cache)

        stages_definition = create_stages(template_standard_uc, c_sys5_uc, c_sys5_ed)
        # Uses IntraStage but the cache is defined in the wrong stage
        sequence_bad_cache = SimulationSequence(
            step_resolution = Hour(24),
            order = Dict(1 => "UC", 2 => "ED"),
            feedforward_chronologies = Dict(("UC" => "ED") => Synchronize(periods = 24)),
            horizons = Dict("UC" => 24, "ED" => 12),
            intervals = Dict(
                "UC" => (Hour(24), Consecutive()),
                "ED" => (Hour(1), Consecutive()),
            ),
            feedforward = Dict(
                ("ED", :devices, :Generators) => SemiContinuousFF(
                    binary_source_stage = PSI.ON,
                    affected_variables = [PSI.ACTIVE_POWER],
                ),
            ),
            cache = Dict(("ED",) => TimeStatusChange(PSY.ThermalStandard, PSI.ON)),
            ini_cond_chronology = IntraStageChronology(),
        )

        sim = Simulation(
            name = "test",
            steps = 1,
            stages = stages_definition,
            stages_sequence = sequence_bad_cache,
            simulation_folder = file_path,
        )
        @test_throws IS.InvalidValue build!(sim)

    end

    @testset "Create Stages with kwargs and custom models" begin
        my_model = JuMP.Model()
        my_model.ext[:PSI_Testing] = 1
        stages_definition_kwargs = Dict(
            "UC" => Stage(
                GenericOpProblem,
                template_basic_uc,
                c_sys5_uc,
                GLPK_optimizer,
                my_model,
            ),
            "ED" => Stage(
                GenericOpProblem,
                template_ed_ptdf,
                c_sys5_ed,
                GLPK_optimizer;
                PTDF = build_PTDF5(),
            ),
        )

        sequence = SimulationSequence(
            step_resolution = Hour(24),
            order = Dict(1 => "UC", 2 => "ED"),
            feedforward_chronologies = Dict(("UC" => "ED") => Synchronize(periods = 24)),
            horizons = Dict("UC" => 24, "ED" => 12),
            intervals = Dict(
                "UC" => (Hour(24), Consecutive()),
                "ED" => (Hour(1), Consecutive()),
            ),
            feedforward = Dict(
                ("ED", :devices, :Generators) => SemiContinuousFF(
                    binary_source_stage = PSI.ON,
                    affected_variables = [PSI.ACTIVE_POWER],
                ),
            ),
            ini_cond_chronology = InterStageChronology(),
        )
        sim = Simulation(
            name = "test",
            steps = 1,
            stages = stages_definition_kwargs,
            stages_sequence = sequence,
            simulation_folder = file_path,
        )
        @test haskey(sim.stages["UC"].internal.psi_container.JuMPmodel.ext, :PSI_Testing)
        @test !isnothing(sim.stages["ED"].internal.psi_container.settings.PTDF)
    end
    @testset "Create Simulation using SOS-PWL cost function" begin
        c_sys5_uc = build_system("c_sys5_pwl_uc")
        c_sys5_ed = build_system("c_sys5_pwl_ed")
        stages_definition_kwargs = Dict(
            "UC" => Stage(
                GenericOpProblem,
                template_pwl_standard_uc,
                c_sys5_uc,
                Cbc_optimizer,
            ),
            "ED" => Stage(GenericOpProblem, template_pwl_ed, c_sys5_ed, Cbc_optimizer),
        )
        sequence = SimulationSequence(
            step_resolution = Hour(24),
            order = Dict(1 => "UC", 2 => "ED"),
            feedforward_chronologies = Dict(("UC" => "ED") => Synchronize(periods = 24)),
            horizons = Dict("UC" => 24, "ED" => 12),
            intervals = Dict(
                "UC" => (Hour(24), Consecutive()),
                "ED" => (Hour(1), Consecutive()),
            ),
            feedforward = Dict(
                ("ED", :devices, :Generators) => SemiContinuousFF(
                    binary_source_stage = PSI.ON,
                    affected_variables = [PSI.ACTIVE_POWER],
                ),
            ),
            ini_cond_chronology = InterStageChronology(),
        )
        sim = Simulation(
            name = "test",
            steps = 2,
            stages = stages_definition_kwargs,
            stages_sequence = sequence,
            simulation_folder = file_path,
        )
        build!(sim)
    end
end

@testset "Test sequence build" begin
    path = joinpath(pwd(), "test_sequence_build")
    !isdir(path) && mkdir(path)
    try
        test_sequence_build(path)
    finally
        @info("removing test files")
        rm(path, force = true, recursive = true)
    end
end

@testset "Test simulation run directory output" begin
    base_path = mktempdir()
    @test PSI._get_output_dir_name(base_path, nothing) == "1"
    mkdir(joinpath(base_path, "1"))
    @test PSI._get_output_dir_name(base_path, nothing) == "2"
    mkdir(joinpath(base_path, "5"))
    @test PSI._get_output_dir_name(base_path, nothing) == "6"
    mkdir(joinpath(base_path, "10"))
    @test PSI._get_output_dir_name(base_path, nothing) == "11"
    @test PSI._get_output_dir_name(base_path, "test") == "test"
end


text:
struct MockMLJ <: Convention end

@testset "void types" begin
    @test scitype(nothing, MockMLJ()) == Nothing
    @test scitype(missing, MockMLJ()) == Missing
    @test scitype([nothing, nothing], MockMLJ()) == AbstractVector{Nothing}
end

@testset "scitype" begin
    X = [1, 2, 3]
    @test scitype(X, MockMLJ()) == AbstractVector{Unknown}

    @test scitype(missing, MockMLJ()) == Missing
    @test scitype((5, 2), MockMLJ()) == Tuple{Unknown,Unknown}
    anyv = Any[5]
    @test scitype(anyv[1], MockMLJ()) == Unknown

    X = [missing, 1, 2, 3]
    @test scitype(X, MockMLJ()) == AbstractVector{Union{Missing, Unknown}}

    Xnm = X[2:end]
    @test scitype(Xnm, MockMLJ()) == AbstractVector{Union{Missing, Unknown}}

    Xm = Any[missing, missing]
    @test scitype(Xm, MockMLJ()) == AbstractVector{Missing}

    @test scitype([missing, missing], MockMLJ()) == AbstractVector{Missing}
end

@testset "scitype2" begin
    ScientificTypesBase.Scitype(::Type{<:Integer}, ::MockMLJ) = Count
    X = [1, 2, 3]
    @test scitype(X, MockMLJ()) == AbstractVector{Count}
    Xm = [missing, 1, 2, 3]
    @test scitype(Xm, MockMLJ()) == AbstractVector{Union{Missing,Count}}
    Xnm = Xm[2:end]
    @test scitype(Xnm, MockMLJ()) == AbstractVector{Union{Missing,Count}}
end

@testset "temporal types" begin
    @test ScientificDate <: ScientificTimeType
    @test ScientificDateTime <: ScientificTimeType
    @test ScientificTime <: ScientificTimeType
end

@testset "compositional" begin
    @test Compositional{3} <: Known
end

@testset "Empty array" begin
    ScientificTypesBase.Scitype(::Type{<:Integer}, ::MockMLJ) = Count
    ScientificTypesBase.Scitype(::Type{Missing}, ::MockMLJ) = Missing
    @test scitype(Int[], MockMLJ()) == AbstractVector{Count}
    @test scitype(Any[], MockMLJ()) == AbstractVector{Unknown}
    @test scitype(Missing[], MockMLJ()) == AbstractVector{Missing}
    @test scitype(Vector{Union{Int,Missing}}(), MockMLJ()) == AbstractVector{Union{Missing,Count}}
end


text:
import Test
import HTTP
import App

function router()
    config = App.Config.getconfig()
    state = App.State.SharedState(config)
    router = App.Router.HttpRouter(state)
    return router
end

@testset "HealthHandler" begin
    @test begin
        r = router()

        req = HTTP.Request("GET", "/health")
        res = HTTP.handle(r, req)

        res.status == 200
    end
end

max_stars_repo_name: DevOpsDemoTF/DevOpsDemo-template-Julia

max_issues_repo_name: DevOpsDemoTF/DevOpsDemo-template-Julia

max_forks_repo_name: DevOpsDemoTF/DevOpsDemo-template-Julia


text:
"""
    mpi_all_to_all(data, comm = mpi_comm(), root = mpi_root)

Declare `data` as mpi-distributed and set communication strategy to all-to-all.

Sets up the [`MPIData`](@ref) structure with [`MPIAllToAll`](@ref) strategy.
"""
function mpi_all_to_all(data, comm = mpi_comm(), root = mpi_root)
    MPI.Initialized() || error("MPI needs to be initialised first.")
    np = MPI.Comm_size(comm)
    id = MPI.Comm_rank(comm)
    s = MPIAllToAll(eltype(storage(data)), np, id, comm)
    return MPIData(data, comm, root, s)
end

"""
     MPIAllToAll

All-to-all communication strategy. The communication works in two steps: first
`MPI.Alltoall!` is used to communicate the number of walkers each rank wants to send to
other ranks, then `MPI.Alltoallv!` is used to send the walkers around.

# Constructor

* `MPIAllToAll(Type{P}, np, id, comm)`: Construct an instance with pair type `P` on
  `np` processes with current rank `id`.

"""
struct MPIAllToAll{P} <: DistributeStrategy
    np::Int32
    id::Int32
    comm::MPI.Comm

    targets::Vector{UInt}
    lenbuffer::MPI.UBuffer{Vector{Cint}} # for communicating numbers of elements
    sendbuffer::MPI.VBuffer{Vector{P}}   # for sending chunks
    recvbuffer::MPI.VBuffer{Vector{P}}   # for receiving chunks
end

function MPIAllToAll(::Type{Pair{K,V}}, np, id, comm) where {K,V}
    P = Pair{K,V}
    datatype = MPI.Datatype(P)
    lenbuf = MPI.UBuffer(zeros(Cint, np), 1, np, MPI.Datatype(Cint))
    sendbuf = MPI.VBuffer(P[], zeros(Cint, np), zeros(Cint, np), datatype)
    recvbuf = MPI.VBuffer(P[], zeros(Cint, np), zeros(Cint, np), datatype)

    return MPIAllToAll{P}(np, id, comm, UInt[], lenbuf, sendbuf, recvbuf)
end

"""
    @swap! arr i j

Swap the `i`-th and `j`-th indices in `arr`.
"""
macro swap!(arr, i, j)
    return quote
        tmp = $(esc(arr))[$(esc(j))]
        $(esc(arr))[$(esc(j))] = $(esc(arr))[$(esc(i))]
        $(esc(arr))[$(esc(i))] = tmp
        $(esc(arr))
    end
end

"""
    sort_and_count!(counts, displs, vec, order, (lo, hi), i_start, j_start)

Sort new spawns by target rank. While this is done, also count the values and calculate the
offsets needed for `MPI.Alltoallv!`.

`counts`, `displs`, `vec`, and `order` are modified in-place. `order` should contain the
values you want to sort vec by. (i.e. `targetrank.(vec, s.np)`)

    sort_and_count!(s::MPIAllToAll)

As above, but operating on the internal buffers of `s`. Note that `s.targets` is expected to
contain the correct values to sort by.
"""
function sort_and_count!(
    counts, displs, vec, ord, (lo, hi), i_start=firstindex(ord), j_start=lastindex(ord)
)
    if i_start == j_start
        # Only one value left. We can stop recursing.
        counts[ord[i_start] + 1] = 1
        displs[ord[i_start] + 1] = i_start - 1
    elseif i_start < j_start
        if lo == hi
            # At this point, the recursion is finished and all values have the same
            # target rank.
            counts[lo + 1] = j_start - i_start + 1
            displs[lo + 1] = i_start - 1
        else
            # Idea:
            # Move from left (i) and right (j). Swap values so that the left part
            # will contain the values that are ≤ mid, while the right part will contain
            # those that are greater.
            # This is essentially equivalent to quicksort, but in this case, we know how
            # to pick a good pivot.
            mid = fld(lo + hi, 2)
            i = i_start
            j = j_start
            @inbounds while i < j
                if ord[i] ≤ mid
                    i += 1
                elseif ord[j] > mid
                    j -= 1
                else
                    @swap! ord i j
                    @swap! vec i j
                end
            end
            # This correction is needed if all values ≤ mid
            if ord[j] ≤ mid
                i += 1
            end
            # Recursively sort subarrays.
            sort_and_count!(counts, displs, vec, ord, (mid + 1, hi), j, j_start)
            sort_and_count!(counts, displs, vec, ord, (lo, mid), i_start, i - 1)
        end
    end
    return nothing
end

function sort_and_count!(s::MPIAllToAll)
    sb = s.sendbuffer
    sb.counts .= 0
    sb.displs .= 0
    sort_and_count!(sb.counts, sb.displs, sb.data, s.targets, (0, s.np - 1))
end

function prepare_send!(s::MPIAllToAll, source)
    buffer = s.sendbuffer.data
    targets = s.targets

    # Copy pairs to send buffer and sort by target rank
    len = length(source)
    if len > 0
        resize!(buffer, len)
        resize!(targets, len)
        @inbounds for (i, p) in enumerate(pairs(source))
            buffer[i] = p
            targets[i] = targetrank(p, s.np)
        end
        sort_and_count!(s)
    end
    copyto!(s.lenbuffer.data, s.sendbuffer.counts)
    return s
end

function prepare_recv!(s::MPIAllToAll)
    counts = s.lenbuffer.data
    resize!(s.recvbuffer.data, sum(counts))
    copyto!(s.recvbuffer.counts, counts)

    displ = 0
    for (r, c) in enumerate(counts)
        s.recvbuffer.displs[r] = displ
        displ += c
    end
    return s
end

function mpi_combine_walkers!(target, source, s::MPIAllToAll)
    prepare_send!(s, source)
    MPI.Barrier(s.comm)
    MPI.Alltoall!(MPI.IN_PLACE, s.lenbuffer, s.comm)
    prepare_recv!(s)
    MPI.Alltoallv!(s.sendbuffer, s.recvbuffer, s.comm)

    for (k, v) in s.recvbuffer.data
        deposit!(target, k, v, nothing)
    end
    return target
end


text:
"""
Check the accuracy between predicted values (Ŷ) and the true value (Y).
"""
function assess_accuracy(Ŷ, Y)
    @assert size(Ŷ)[1] == size(Y)[1]
    return Ŷ[findmax(Y)[2]]
end # function

max_stars_repo_name: Squalm/NeuralNetworks

max_issues_repo_name: Squalm/NeuralNetworks

max_forks_repo_name: Squalm/NeuralNetworks


text:
using Test
include("../../src/contrib/taxeffect.jl")

taxes = rand(10) * 100

m = getuniformtaxmodel()
update_param!(m, :uniformtax, taxes)
run(m)

for ii in 1:10
    @test all(m[:AbatementCostsCO2, :mc_marginalcost][ii, :] .≈ taxes[ii])
    @test all(m[:AbatementCostsCH4, :mc_marginalcost][ii, :] .≈ taxes[ii])
    @test all(m[:AbatementCostsN2O, :mc_marginalcost][ii, :] .≈ taxes[ii])
    @test all(m[:AbatementCostsLin, :mc_marginalcost][ii, :] .≈ taxes[ii])
end