text:
using VennEuler, Test, Random

# just use @test macro
# run with:
# reload("src/VennEuler.jl"); include("test/test.jl")

# set up an initial test of 3 circles with random set membership

# data is a boolean matrix of N rows/observations, with M columns/sets,
# and any number of set memberships for each observation (although often one
# or more in practice).
randdata = [
  true   true  false
  true  false   true
  true  false  false
  true  false  false
  true  false  false
 false  false   true
  true  false   true
  true  false   true
 false  false   true
  true  false   true
 false  false  false
 false   true  false
 false  false  false
  true  false  false
  true  false   true
 false  false   true
 false   true   true
 false   true  false
  true  false   true
  true   true  false ]

setlabels = ["A", "B", "C"]
# specification is either a single object (if all sets are treated the same)
# or a list of objects (to treat them differently)
spec = EulerSpec(:circle, [NaN, NaN], [0, 0])
specdefault = EulerSpec()
@test isequal(spec.shape, specdefault.shape)
@test isequal(spec.clamp, specdefault.clamp)
@test_throws ErrorException EulerSpec(:tesseract)
@test_throws AssertionError EulerSpec(:circle, [1, 1, 1], [0, 0, 0])

spec2 = EulerSpec(:circle, [.5, NaN], [0, 0])
specs1 = [deepcopy(spec), deepcopy(spec), spec2]
VennEuler.update_statepos!(specs1)
@test isequal(specs1[1].statepos, [1, 2])
@test isequal(specs1[2].statepos, [3, 4])

ss, bb = VennEuler.compute_shape_sizes(specs1, vec(sum(randdata, dims=1)), .5)
@test ss ≈ [0.265962,0.171677,0.242789] atol=.001

# make sure DisjointSets can be constructed
ds1 = VennEuler.DisjointSet(randdata, setlabels)
@test ds1.counts == [2,3,2,1,4,6,2,0]

# make an Euler object
eo = make_euler_object(setlabels, randdata, spec, sizesum=.5) # test shortcut
@test eo.lb ≈ [0.265962, 0.265962, 0.171677, 0.171677, 0.242789, 0.242789] atol=.001
es = random_state(eo)
@test all(eo.lb .<= es .<= eo.ub)

@test VennEuler.get_params(es, specs1[1]) == tuple(es[1:2]...)

# tests for bitmap operations
bmc1 = VennEuler.make_bitmap_circle(.5, .5, .2, 20)
bmc2 = VennEuler.make_bitmap_circle(.5, .5, .001, 20)
@test sum(bmc1) == 57
println("")
VennEuler.showbitmap(bmc1)
@test sum(bmc2) == 1

spec_sq = EulerSpec(:square, [NaN, NaN], [0, 0])
bmsq1 = VennEuler.make_bitmap([.5, .5], .2, spec_sq, 20)
@test sum(bmsq1) == 81
println("")
VennEuler.showbitmap(bmsq1)

spec_tr = EulerSpec(:triangle, [NaN, NaN], [0, 0])
bmtr1 = VennEuler.make_bitmap([.5, .5], .15, spec_tr, 20)
#@show bmtr1 #sum(bmsq1) == 81
println("")
VennEuler.showbitmap(bmtr1)

spec_rect = EulerSpec(:rectangle, [NaN, NaN, NaN], [0, 0, 0])
bmrect1 = VennEuler.make_bitmap([.5, .5, 1], .2, spec_rect, 20)
bmrect2 = VennEuler.make_bitmap([.5, .5, 0], .2, spec_rect, 20)
println("")
VennEuler.showbitmap(bmrect1)
println("")
VennEuler.showbitmap(bmrect2)

println("")

# evaluating Euler states
simplelabels = ["A", "B"]
simpledata = Bool[1 0; 1 1; 0 1]

eo2 = make_euler_object(simplelabels, simpledata, spec, sizesum=.5)
@show VennEuler.eval_euler_state(eo2, [.33, .5, .66, .5])
@show eo2.evalfn([.33, .5, .66, .5], [])
@time VennEuler.eval_euler_state(eo2, [.38, .5, .62, .5])

# optimization
(minf,minx,ret) = optimize(eo2, random_state(eo2))
println("got $minf at $minx (returned $ret)")
#@test ret == :FTOL_REACHED

eo = make_euler_object(setlabels, randdata, [spec, spec_tr, spec_sq], sizesum=.5)
(minf,minx,ret) = optimize(eo, random_state(eo), ftol=1/10000)
println("got $minf at $minx (returned $ret)")
#@test ret == :FTOL_REACHED

spec_fixed = EulerSpec(:circle, [0.5, 0.5], [0, 0])
eo = make_euler_object(setlabels, randdata, [spec_fixed, spec_tr, spec_rect], sizesum=.2)
# TODO: better feedback if sizesum is too big
(minf,minx,ret) = optimize(eo, random_state(eo), ftol=1/10000)
println("got $minf at $minx (returned $ret)")
@test isequal(minx[1:2], spec_fixed.clamp)

(minf,minx,ret) = optimize_iteratively(eo, random_state(eo), ftol=1/10000, verbose=true)
println("got $minf at $minx (returned $ret)")

# output
render("test1.svg", eo, minx, verbose=1)


text:
using Revise, ProgressMeter, OffsetArrays, Kinetic
#include("/home/vavrines/Coding/SKS.jl/src/SKS.jl")
#using .SKS
begin
    using Reexport
    using OffsetArrays
    using Plots
    using JLD2
    using PolyChaos

    include("/home/vavrines/Coding/Langevin.jl/src/uq.jl")
    include("/home/vavrines/Coding/Langevin.jl/src/kinetic.jl")
    include("/home/vavrines/Coding/Langevin.jl/src/initialize.jl")
    include("/home/vavrines/Coding/Langevin.jl/src/flux.jl")
    include("/home/vavrines/Coding/Langevin.jl/src/out.jl")
    include("/home/vavrines/Coding/Langevin.jl/src/solver.jl")
    include("/home/vavrines/Coding/Langevin.jl/src/step.jl")
end

begin
    case = "weibel"
    space = "1d3f2v"
    nSpecies = 2
    interpOrder = 2
    limiter = "vanleer"
    cfl = 0.3
    maxTime = 0.1

    # physical space
    x0 = 0
    x1 = 1
    nx = 100
    pMeshType = "uniform"
    nxg = 2

    # velocity space
    vMeshType = "rectangle"
    umin = -5
    umax = 5
    nu = 28
    nug = 0

    vmin = -5
    vmax = 5
    nv = 28
    nvg = 0

    # random space
    uqMethod = "collocation"
    nr = 1
    nRec = 3
    opType = "uniform"
    parameter1 = 0.95
    parameter2 = 1.05

    # gas
    knudsen = 0.000001
    mach = 0.0
    prandtl = 1
    inK = 0

    mi = 1
    ni = 0.5
    #me = 0.001
    me = 0.0005446623
    ne = 0.5
    lD = 0.01
    rL = 0.003

    # electromagnetic field
    sol = 100
    echi = 1
    bnu = 1
end

begin
    γ = heat_capacity_ratio(inK, 3)
    set = Setup(case, space, nSpecies, interpOrder, limiter, cfl, maxTime)
    pSpace = PSpace1D(x0, x1, nx, pMeshType, nxg)

    ue0 = umin * sqrt(mi / me)
    ue1 = umax * sqrt(mi / me)
    ve0 = vmin * sqrt(mi / me)
    ve1 = vmax * sqrt(mi / me)
    kne = knudsen * (me / mi)

    vSpace =
        MVSpace2D(umin, umax, ue0, ue1, nu, vmin, vmax, ve0, ve1, nv, vMeshType, nug, nvg)
    plasma = Plasma2D(
        [knudsen, kne],
        mach,
        prandtl,
        inK,
        γ,
        mi,
        ni,
        me,
        ne,
        lD,
        rL,
        sol,
        echi,
        bnu,
    )

    begin
        # upstream
        primL = zeros(5, 2)
        primL[1, 1] = 1.0 * mi
        primL[2, 1] = 0.0
        primL[3, 1] = 0.0
        primL[4, 1] = 0.0
        primL[5, 1] = mi / 1.0
        primL[1, 2] = 1.0 * me
        primL[2, 2] = 0.0
        primL[3, 2] = 0.0
        primL[4, 2] = 0.0
        primL[5, 2] = me / 1.0

        wL = mixture_prim_conserve(primL, γ)
        h0L = mixture_maxwellian(vSpace.u, vSpace.v, primL)

        h1L = similar(h0L)
        h2L = similar(h0L)
        for j in axes(h0L, 3)
            h1L[:, :, j] .= primL[4, j] .* h0L[:, :, j]
            h2L[:, :, j] .= (primL[4, j]^2 + 1.0 / (2.0 * primL[end, j])) .* h0L[:, :, j]
        end

        EL = zeros(3)
        BL = zeros(3)
        BL[1] = 0.75
        BL[2] = 1.0

        # downstream
        primR = zeros(5, 2)
        primR[1, 1] = 0.125 * mi
        primR[2, 1] = 0.0
        primR[3, 1] = 0.0
        primR[4, 1] = 0.0
        primR[5, 1] = mi * 1.25
        primR[1, 2] = 0.125 * me
        primR[2, 2] = 0.0
        primR[3, 2] = 0.0
        primR[4, 2] = 0.0
        primR[5, 2] = me * 1.25

        wR = mixture_prim_conserve(primR, γ)
        h0R = mixture_maxwellian(vSpace.u, vSpace.v, primR)

        h1R = similar(h0R)
        h2R = similar(h0R)
        for j in axes(h0R, 3)
            h1R[:, :, j] .= primR[4, j] .* h0R[:, :, j]
            h2R[:, :, j] .= (primR[4, j]^2 + 1.0 / (2.0 * primR[end, j])) .* h0R[:, :, j]
        end

        ER = zeros(3)
        BR = zeros(3)
        BR[1] = 0.75
        BR[2] = -1.0

        lorenzL = zeros(3, 2)
        lorenzR = zeros(3, 2)
        bcL = zeros(5, 2)
        bcR = zeros(5, 2)
    end

    ib = IB3F(
        wL,
        primL,
        h0L,
        h1L,
        h2L,
        bcL,
        EL,
        BL,
        lorenzL,
        wR,
        primR,
        h0R,
        h1R,
        h2R,
        bcR,
        ER,
        BR,
        lorenzR,
    )

    outputFolder = pwd()

    ks = SolverSet(set, pSpace, vSpace, plasma, ib, outputFolder)
    KS = ks

    ctr = OffsetArray{ControlVolume1D3F}(undef, axes(KS.pSpace.x, 1))
    face = Array{Interface1D3F}(undef, KS.pSpace.nx + 1)

    uq = UQ1D(nr, nRec, parameter1, parameter2, opType, uqMethod)
end

begin
    idx0 = (eachindex(pSpace.x)|>collect)[1]
    idx1 = (eachindex(pSpace.x)|>collect)[end]

    # upstream
    primL = zeros(5, uq.op.quad.Nquad, 2)
    for j in axes(primL, 2)
        primL[:, j, :] .= KS.ib.primL
    end

    wL = uq_prim_conserve(primL, KS.gas.γ, uq)
    h0L, h1L, h2L = uq_maxwellian(KS.vSpace.u, KS.vSpace.v, primL, uq)

    EL = zeros(3, uq.op.quad.Nquad)
    BL = zeros(3, uq.op.quad.Nquad)
    for j in axes(BL, 2)
        BL[:, j] .= KS.ib.BL
    end

    # downstream
    primR = zeros(5, uq.op.quad.Nquad, 2)
    for j in axes(primL, 2)
        primR[:, j, :] .= KS.ib.primR
    end

    wR = uq_prim_conserve(primR, KS.gas.γ, uq)
    h0R, h1R, h2R = uq_maxwellian(KS.vSpace.u, KS.vSpace.v, primR, uq)

    ER = zeros(3, uq.op.quad.Nquad)
    BR = zeros(3, uq.op.quad.Nquad)
    for j in axes(BR, 2)
        BR[:, j] .= KS.ib.BR
    end

    lorenz = zeros(3, uq.op.quad.Nquad, 2)

    for i in eachindex(ctr)
        if i <= KS.pSpace.nx ÷ 2
            ctr[i] = ControlVolume1D3F(
                KS.pSpace.x[i],
                KS.pSpace.dx[i],
                wL,
                primL,
                h0L,
                h1L,
                h2L,
                EL,
                BL,
                lorenz,
            )
        else
            ctr[i] = ControlVolume1D3F(
                KS.pSpace.x[i],
                KS.pSpace.dx[i],
                wR,
                primR,
                h0R,
                h1R,
                h2R,
                ER,
                BR,
                lorenz,
            )
        end
    end

    face = Array{Interface1D3F}(undef, KS.pSpace.nx + 1)
    for i = 1:KS.pSpace.nx+1
        face[i] = Interface1D3F(wL, h0L, EL)
    end
end

iter = 0
res = zeros(5, 2)
simTime = 0.0
dt = timestep(KS, ctr, simTime, uq)
#dt = SKS.timestep(KS, ctr, simTime, uq)
nt = Int(floor(ks.set.maxTime / dt))

res = zeros(5, 2)
@showprogress for iter = 1:nt
    #dt = timestep(KS, ctr, simTime, uq)
    #    KitBase.reconstruct!(KS, ctr)

    #evolve!(KS, uq, ctr, face, dt)

    @inbounds Threads.@threads for i in eachindex(face)
        uqflux_flow!(KS, ctr[i-1], face[i], ctr[i], dt, mode = :kfvs)

        for j = 1:uq.op.quad.Nquad
            femL = @view face[i].femL[:, j]
            femR = @view face[i].femR[:, j]

            flux_em!(
                femL,
                femR,
                ctr[i-2].E[:, j],
                ctr[i-2].B[:, j],
                ctr[i-1].E[:, j],
                ctr[i-1].B[:, j],
                ctr[i].E[:, j],
                ctr[i].B[:, j],
                ctr[i+1].E[:, j],
                ctr[i+1].B[:, j],
                ctr[i-1].ϕ[j],
                ctr[i].ϕ[j],
                ctr[i-1].ψ[j],
                ctr[i].ψ[j],
                ctr[i-1].dx,
                ctr[i].dx,
                KS.gas.A1p,
                KS.gas.A1n,
                KS.gas.D1,
                KS.gas.sol,
                KS.gas.χ,
                KS.gas.ν,
                dt,
            )
        end

    end

    update!(KS, uq, ctr, face, dt, res)

    #iter += 1
    #t += dt

    #if iter%1000 == 0
    #    println("iter: $(iter), time: $(simTime), dt: $(dt), res: $(res[1:5,1:2])")
    #end

    #if simTime > KS.set.maxTime || maximum(res) < 1.e-7
    #    break
    #end
end

sol = zeros(ks.pSpace.nx, 10, 2)
for i = 1:ks.pSpace.nx
    sol[i, 1, 1] = ctr[i].prim[1, 1, 1]
    sol[i, 1, 2] = ctr[i].prim[1, 1, 2] / ks.gas.me
    sol[i, 2:4, 1] .= ctr[i].prim[2:4, 1, 1]
    sol[i, 2:4, 2] .= ctr[i].prim[2:4, 1, 2]
    sol[i, 5, 1] = 1.0 / ctr[i].prim[5, 1, 1]
    sol[i, 5, 2] = ks.gas.me / ctr[i].prim[5, 1, 2]

    sol[i, 6, 1] = ctr[i].B[2, 1]
    sol[i, 6, 2] = ctr[i].E[1, 1]
end
using Plots
plot(ks.pSpace.x[1:ks.pSpace.nx], sol[:, 1, 1])
plot!(ks.pSpace.x[1:ks.pSpace.nx], sol[:, 1, 2])

plot(ks.pSpace.x[1:ks.pSpace.nx], sol[:, 6, 1])
plot!(ks.pSpace.x[1:ks.pSpace.nx], sol[:, 6, 2])

@save "sol.jld2" ks uq ctr

ext: jl


text:
struct FixedEffectSolverLSMR{T} <: AbstractFixedEffectSolver{T}
	m::FixedEffectLinearMap{T}
	b::AbstractVector{T}
	r::AbstractVector{T}
	x::FixedEffectCoefficients{T}
	v::FixedEffectCoefficients{T}
	h::FixedEffectCoefficients{T}
	hbar::FixedEffectCoefficients{T}
end

function AbstractFixedEffectSolver{T}(fes::Vector{<:FixedEffect}, weights::AbstractWeights, ::Type{Val{:lsmr}}) where {T}
	m = FixedEffectLinearMap{T}(fes, weights, Val{:lsmr})
	b = zeros(T, length(weights))
	r = zeros(T, length(weights))
	x = FixedEffectCoefficients([zeros(T, fe.n) for fe in fes])
	v = FixedEffectCoefficients([zeros(T, fe.n) for fe in fes])
	h = FixedEffectCoefficients([zeros(T, fe.n) for fe in fes])
	hbar = FixedEffectCoefficients([zeros(T, fe.n) for fe in fes])
	return FixedEffectSolverLSMR(m, b, r, x, v, h, hbar)
end

function solve_residuals!(r::AbstractVector, feM::FixedEffectSolverLSMR{T}; tol::Real = sqrt(eps(T)), maxiter::Integer = 100_000) where {T}
	copyto!(feM.r, r)
	feM.r .*=  feM.m.sqrtw
	fill!(feM.x, 0)
	copyto!(feM.b, feM.r)
	x, ch = lsmr!(feM.x, feM.m, feM.b, feM.v, feM.h, feM.hbar; atol = tol, btol = tol, maxiter = maxiter)
	mul!(feM.r, feM.m, feM.x, -1.0, 1.0)
	feM.r ./=  feM.m.sqrtw
	copyto!(r, feM.r)
	return r, div(ch.mvps, 2), ch.isconverged
end

function solve_residuals!(X::AbstractMatrix, feM::FixedEffectSolverLSMR; kwargs...)
	iterations = Int[]
	convergeds = Bool[]
	for j in 1:size(X, 2)
		_, iteration, converged = solve_residuals!(view(X, :, j), feM; kwargs...)
		push!(iterations, iteration)
		push!(convergeds, converged)
	end
	return X, iterations, convergeds
end

function solve_coefficients!(r::AbstractVector, feM::FixedEffectSolverLSMR{T}; tol::Real = sqrt(eps(T)), maxiter::Integer = 100_000) where {T}
	copyto!(feM.b, r)
	feM.b .*=  feM.m.sqrtw
	fill!(feM.x, 0)
	x, ch = lsmr!(feM.x, feM.m, feM.b, feM.v, feM.h, feM.hbar; atol = tol, btol = tol, maxiter = maxiter)
	for (x, scale) in zip(feM.x.x, feM.m.colnorm)
		x ./=  scale
	end
	x = Vector{eltype(r)}[x for x in feM.x.x]
	full(normalize!(x, feM.m.fes; tol = tol, maxiter = maxiter), feM.m.fes), div(ch.mvps, 2), ch.isconverged
end


text:
scalar(a::Number) = a
scalar(A) = A[1]

## Taken from version 0.4
cld{T<:Integer }(x::T, y::T) = div(x,y)+(!signbit(x$y)&(rem(x,y)!=0))


text:
# Autogenerated wrapper script for libidn2_jll for armv7l-linux-gnueabihf
export idn2, libidn2

using Libiconv_jll
using libunistring_jll
JLLWrappers.@generate_wrapper_header("libidn2")
JLLWrappers.@declare_executable_product(idn2)
JLLWrappers.@declare_library_product(libidn2, "libidn2.so.0")
function __init__()
    JLLWrappers.@generate_init_header(Libiconv_jll, libunistring_jll)
    JLLWrappers.@init_executable_product(
        idn2,
        "bin/idn2",
    )

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

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
initvalues(model::GMHPhotoReceptor.PhotoReceptorModel) = GeneralizedMetropolisHastings.initvalues(trait(:initialize,:prior),parameters(model),Float64)

#Iterate evaluate n times and store all currents + all logposteriors
function evaluate1(model::GMHPhotoReceptor.PhotoReceptorModel,niterations::Integer;paramvals=nothing)
    if paramvals == nothing
        paramvals = initvalues(model)
    end
    ntimesteps = GeneralizedMetropolisHastings.numvalues(model.photons)
    macrocurrents = zeros(Float32,ntimesteps,niterations)
    logposteriors = zeros(Float32,niterations)
    lprior = logprior(parameters(model),paramvals,Float32)
    for i=1:niterations
        macrocurrent =  evaluate!(model,paramvals)
        copy!(macrocurrents,(i-1)*ntimesteps+1,macrocurrent,1,ntimesteps)
        logposteriors[i] = loglikelihood(model,macrocurrent) + lprior
    end
    paramvals,macrocurrents,logposteriors
end

#Iterate evaluate n times and store mean and std current + all logposteriors
function evaluate2(model::GMHPhotoReceptor.PhotoReceptorModel,niterations::Integer;paramvals=nothing,macrocurrentstats=nothing)
    ntimesteps = GeneralizedMetropolisHastings.numvalues(model.photons)
    if paramvals == nothing
        paramvals = initvalues(model)
    end
    if macrocurrentstats == nothing
        macrocurrentstats = StreamingStats(ntimesteps,Float32)
    end
    logposteriors= zeros(Float32,niterations)
    lprior = logprior(parameters(model),paramvals,Float32)
    for i=1:niterations
        macrocurrent =  evaluate!(model,paramvals)
        update!(macrocurrentstats,macrocurrent)
        logposteriors[i] = loglikelihood(model,macrocurrent) + lprior
    end
    paramvals,macrocurrentstats,logposteriors
end

#Run evaluate2 for 4 different parameter values and return all mean and std of macrocurrents + logposteriors
function evaluate3(model::GMHPhotoReceptor.PhotoReceptorModel,niterations::Integer,nparamsets::Integer)
    ntimesteps = GeneralizedMetropolisHastings.numvalues(model.photons)
    paramvals = zeros(GeneralizedMetropolisHastings.numparas(model),nparamsets)
    meancurrents = zeros(Float32,ntimesteps,nparamsets)
    stdcurrents = zeros(meancurrents)
    logposteriors = zeros(Float32,niterations,nparamsets)
    r = pmap((n)->evaluate2(model,niterations),1:nparamsets)
    for p=1:nparamsets
        paramvals[:,p] = r[p][1]
        meancurrents[:,p] = mean(r[p][2])
        stdcurrents[:,p] = std(r[p][2])
        logposteriors[:,p] = r[p][3]
    end
    paramvals,meancurrents,stdcurrents,logposteriors
end

function gradient1(model::GMHPhotoReceptor.PhotoReceptorModel,niterations::Integer,gradienteps::Number=1e-4;paramvals=nothing)
    ntimesteps = GeneralizedMetropolisHastings.numvalues(model.photons)
    nparamsets = GeneralizedMetropolisHastings.numparas(model) + 1
    gradparamvals = repmat(initvalues(model),1,nparamsets)
    for i=2:nparamsets
        gradparamvals[i-1,i]+=gradienteps
    end
    meancurrents = zeros(Float32,ntimesteps,nparamsets)
    stdcurrents = zeros(meancurrents)
    logposteriors = zeros(Float32,niterations,nparamsets)
    r = pmap((n)->evaluate2(model,niterations;paramvals=gradparamvals[:,n]),1:nparamsets)
    for p=1:nparamsets
        meancurrents[:,p] = mean(r[p][2])
        stdcurrents[:,p] = std(r[p][2])
        logposteriors[:,p] = r[p][3]
    end
    gradparamvals,meancurrents,stdcurrents,logposteriors
end


text:
using OpenMDAO
using PyCall
using CCBlade
using PyPlot

om = pyimport("openmdao.api")
julia_comps = pyimport("omjl.julia_comps")
pyccblade = pyimport("ccblade.ccblade_jl")
pyccblade_geom = pyimport("ccblade.geometry")

num_nodes = 1
num_blades = 3
num_radial = 15

af_filename = "airfoils/mh117.dat"
af = af_from_file(af_filename, use_interpolations_jl=true)

num_cp = 6
chord = 10.0
theta = collect(LinRange(65., 25., num_cp)) .* pi/180.0
pitch = 0.0

hub_diameter = 30.  # cm
prop_diameter = 150.  # cm
c0 = sqrt(1.4*287.058*300.0)  # meters/second
rho0 = 1.4*98600.0/(c0*c0)  # kg/m^3
omega = 236.0

prob = om.Problem()

comp = om.IndepVarComp()
comp.add_output("rho", val=rho0, shape=num_nodes, units="kg/m**3")
comp.add_output("mu", val=1., shape=num_nodes, units="N/m**2*s")
comp.add_output("asound", val=c0, shape=num_nodes, units="m/s")
comp.add_output("v", val=77.2, shape=num_nodes, units="m/s")
comp.add_output("alpha", val=0., shape=num_nodes, units="rad")
comp.add_output("incidence", val=0., shape=num_nodes, units="rad")
comp.add_output("precone", val=0., units="deg")
comp.add_output("omega", val=omega, shape=num_nodes, units="rad/s")
comp.add_output("hub_diameter", val=hub_diameter, shape=num_nodes, units="cm")
comp.add_output("prop_diameter", val=prop_diameter, shape=num_nodes, units="cm")
comp.add_output("pitch", val=pitch, shape=num_nodes, units="rad")
comp.add_output("chord_dv", val=chord, shape=num_cp, units="cm")
comp.add_output("theta_dv", val=theta, shape=num_cp, units="rad")
prob.model.add_subsystem("indep_var_comp", comp, promotes=["*"])

comp = pyccblade_geom.GeometryGroup(num_nodes=num_nodes, num_cp=num_cp, num_radial=num_radial)
prob.model.add_subsystem(
    "geometry_group", comp,
    promotes_inputs=["hub_diameter", "prop_diameter", "chord_dv",
                     "theta_dv", "pitch"],
    promotes_outputs=["radii", "dradii", "chord", "theta"])

group = om.Group()

comp = om.ExecComp("hub_radius = 0.5*hub_diameter",
                   hub_radius=Dict("value" => 0.1, "units" => "m"),
                   hub_diameter=Dict("units" => "m"))
group.add_subsystem("hub_radius_comp", comp, promotes=["*"])

comp = om.ExecComp("prop_radius = 0.5*prop_diameter",
                   prop_radius=Dict("value" => 1.0, "units" => "m"),
                   prop_diameter=Dict("units" => "m"))
group.add_subsystem("prop_radius_comp", comp, promotes=["*"])

comp = om.ExecComp("Vx = v*cos(precone)",
                   v=Dict("units" => "m/s",
                          "shape" => (num_nodes, num_radial),
                          "src_indices" => repeat(0:num_nodes-1, num_nodes, num_radial)),
                   precone=Dict("units" => "rad"),
                   Vx=Dict("units" => "m/s", "shape" => (num_nodes, num_radial)))
group.add_subsystem("Vx_comp", comp, promotes=["*"])

comp = om.ExecComp("Vy = omega*radii*cos(precone)",
                   omega=Dict("units" => "rad/s",
                          "shape" => (num_nodes, 1),
                          "flat_src_indices" => collect(0:num_nodes-1)),
                   radii=Dict("units" => "m", "shape" => (1, num_radial)),
                   precone=Dict("units" => "rad"),
                   Vy=Dict("units" => "m/s", "shape" => (num_nodes, num_radial)))
group.add_subsystem("Vy_comp", comp, promotes=["*"])

ccblade_residual_comp_data = CCBlade.CCBladeResidualComp(
    num_nodes=num_nodes, num_radial=num_radial, af=af, B=num_blades, turbine=false, debug_print=true)
comp = julia_comps.JuliaImplicitComp(julia_comp_data=ccblade_residual_comp_data)
comp.linear_solver = om.DirectSolver(assemble_jac=true)
# comp.nonlinear_solver = om.NewtonSolver(solve_subsystems=true, iprint=2, err_on_non_converge=true)
group.add_subsystem("ccblade_residual_comp", comp,
                    promotes_inputs=[("r", "radii"), "chord", "theta", "Vx",
                                     "Vy", "rho", "mu", "asound", "pitch",
                                     ("Rhub", "hub_radius"), ("Rtip", "prop_radius"),
                                     "precone"],
                    promotes_outputs=["Np", "Tp"])

comp = pyccblade.CCBladeThrustTorqueComp(num_nodes=num_nodes,
                                         num_radial=num_radial, B=num_blades)
group.add_subsystem("ccblade_torquethrust_comp", comp,
                    promotes_inputs=["radii", "dradii", "Np", "Tp"],
                    promotes_outputs=["thrust", "torque"])

comp = om.ExecComp("efficiency = (thrust*v)/(torque*omega)",
                   thrust=Dict("units" => "N", "shape" => num_nodes),
                   v=Dict("units" => "m/s", "shape" => num_nodes),
                   torque=Dict("units" => "N*m", "shape" => num_nodes),
                   omega=Dict("units" => "rad/s", "shape" => num_nodes),
                   efficiency=Dict("shape" => num_nodes))
group.add_subsystem("efficiency_comp", comp,
                    promotes_inputs=["*"],
                    promotes_outputs=["*"])

group.linear_solver = om.DirectSolver(assemble_jac=true)

prob.model.add_subsystem(
    "ccblade_group", group,
    promotes_inputs=["radii", "dradii", "chord", "theta", "rho", "mu",
                     "asound", "v", "precone", "omega", "hub_diameter",
                     "prop_diameter", "pitch"],
    promotes_outputs=["thrust", "torque", "efficiency"])

prob.model.add_design_var("chord_dv", lower=1., upper=20.,
                          scaler=5e-2)
prob.model.add_design_var("theta_dv",
                          lower=20.0*pi/180., upper=90.0*pi/180.0)

prob.model.add_objective("efficiency", scaler=-1.,)
prob.model.add_constraint("thrust", equals=700.0, scaler=1e-3,
                          indices=collect(0:num_nodes-1))
prob.driver = om.pyOptSparseDriver(optimizer="SNOPT")

prob.setup()
prob.final_setup()
prob.run_driver()

node = 1
radii = prob.get_val("radii", units="m")[node, :]
ccblade_normal_load = prob.get_val(
    "ccblade_group.Np", units="N/m")[node, :]*num_blades
ccblade_circum_load = prob.get_val(
    "ccblade_group.Tp", units="N/m")[node, :]*num_blades

fig, ax = plt.subplots()
plot(radii, ccblade_normal_load, label="CCBlade.jl")
xlabel("blade element radius, m")
ylabel("normal load, N/m")
legend()
savefig("normal_load.png")

fig, ax = plt.subplots()
plot(radii, ccblade_circum_load, label="CCBlade.jl")
xlabel("blade element radius, m")
ylabel("circumferential load, N/m")
legend()
savefig("circum_load.png")


text:
module TestDeclarations

using Unitful
using BenchmarkTools

struct Var
	var::NamedTuple
	Var(; kwargs...) = new((;kwargs...) )
end

struct Model
	model::NamedTuple
	Model(; kwargs...) = new((;kwargs...) )
end

struct Map
	map::NamedTuple
	Map(; kwargs...) = new((;kwargs...) )
end

Pin = Model( v = Var(potential=true, nominal=10), i = Var(;flow=true) )

@show Pin

Input(; kwargs...) = Var(;input=true, kwargs...)
Output(; kwargs...) = Var(;output=true, kwargs...)
Potential(; kwargs...) = Var(;potential=true, kwargs...)
Flow(; kwargs...) = Var(;flow=true, kwargs...)

using Base.Meta: isexpr
using OrderedCollections: OrderedDict
using Unitful
using ModiaBase.Symbolic: removeBlock

function showModel(m, level=0)
    println("(")
    level += 1
    for (k, v) in zip(keys(m), m)
        if typeof(v) <: NamedTuple
            print("  "^level, k, " = ")
            showModel(v, level)
        else
            println("  "^level, k, " = ", removeBlock(v), ",")
        end
    end
    println("  "^(level-1), "),")
end

macro showModel(model)
    modelName = string(model)
    mod = :( print($modelName, " = Model"); showModel($model); println() )
    return esc(mod)
end

global logMerge = true

function setLogMerge(val)
    global logMerge
    logMerge = val
end

recursiveMerge(x, ::Nothing) = x
recursiveMerge(x, y) = y
recursiveMerge(x::Expr, y::Expr) = begin dump(x); dump(y); Expr(x.head, x.args..., y.args...) end
recursiveMerge(x::Expr, y::Tuple) = begin x = copy(x); xargs = x.args; xargs[y[2]] = y[3]; Expr(x.head, xargs...) end

function recursiveMerge(nt1::NamedTuple, nt2::NamedTuple)
    all_keys = union(keys(nt1), keys(nt2))
    gen = Base.Generator(all_keys) do key
        v1 = get(nt1, key, nothing)
        v2 = get(nt2, key, nothing)
        key => recursiveMerge(v1, v2)
    end
    return (; gen...)
end

function mergeModels(m1::NamedTuple, m2::NamedTuple, env=Symbol())
    mergedModels = OrderedDict{Symbol,Any}(pairs(m1))
    for (k,v) in collect(pairs(m2))
        if typeof(v) <: NamedTuple
            if k in keys(mergedModels) && ! (:_redeclare in keys(v))
                if logMerge; print("In $k: ") end
                m = mergeModels(mergedModels[k], v, k)
                mergedModels[k] = m
            elseif :_redeclare in keys(v)
                if logMerge; println("Redeclaring: $k = $v") end
                mergedModels[k] = v
            elseif nothing in values(v)  # TODO: Refine

            else
                if !(:_redeclare in keys(mergedModels))
                    if logMerge; println("Adding: $k = $v") end
                end
                mergedModels[k] = v
            end
        elseif v === nothing
            if logMerge; println("Deleting: $k") end
            delete!(mergedModels, k)
        else
            if logMerge
                if k in keys(mergedModels)
                    println("Changing: $k = $(mergedModels[k]) to $k = $v")
                elseif !(:_redeclare in keys(mergedModels))
                    println("Adding: $k = $v")
                end
            end
            mergedModels[k] = v
        end
    end
    return (; mergedModels...) # Transform OrderedDict to named tuple
end

mergeModels(m1::Model, m2::Model, env=Symbol()) = recursiveMerge(m1.model, m2.model)
mergeModels(m1::Model, m2::Map, env=Symbol()) = recursiveMerge(m1.model, m2.map)

#Base.:|(m::Model, n::Model) =  recursiveMerge(m.model, n.model)
#Base.:|(m::Model, n::Map) = recursiveMerge(m.model, n.map)
#Base.:|(m, n::Map) = recursiveMerge(m, n.map)
Base.:|(m::Model, n::Model) =  (:MergeModel, m.model, n.model)
Base.:|(m::Model, n::Map) =  (:MergeMap, m.model, n.map)
Base.:|(m, n::Map) =  (:MergeMap, m, n.map)

Redeclare = ( _redeclare = true, )
Replace(i, e) = (:replace, i, e)
Final(x) = x

# ----------------------------------------------------------

@time Pin = Model( v = Potential(), i = Flow() )

OnePort = Model( p = Pin, n = Pin, equations = :[
        v = p.v - n.v
        0 = p.i + n.i
        i = p.i ]
    )

@showModel(OnePort.model)

@time Resistor = OnePort | Model( R = Final(1.0u"Ω"), equations = Replace(1, :( R*i = v ) ))

@showModel(Resistor)

Capacitor = OnePort | Model( C = 1.0u"F", v = Map(init=0.0u"V"), equations = :[ C*der(v) = i ] )

Inductor = OnePort | Model( L = 1.0u"H", i = Map(init=0.0u"A"), equations = :[ L*der(i) = v ] )

ConstantVoltage = OnePort | Model( V = 1.0u"V", equations = :[ v = V ] )

Ground = Model( p = Pin, equations = :[ p.v = 0.0u"V" ] )

@time Filter = Model(
    R = Resistor | Map(R=0.5u"Ω"),
    C = Capacitor | Map(C=2.0u"F", v=Map(init=0.1u"V")),
    V = ConstantVoltage | Map(V=10.0u"V"),
    ground = Ground,
    connect = :[
      (V.p, R.p)
      (R.n, C.p)
      (C.n, V.n, ground.p)
    ]
)

@showModel(Filter.model)

end


text:
# This converts the ground-truth JSON files to the Julia repr format so
# we can use that without requiring a JSON parser during testing.

using JSON

const testfiles =  joinpath(@__DIR__, "..", "testfiles")

function convert_json_files()
    for folder in ("invalid", "valid")
        for file in readdir(joinpath(testfiles, folder); join=true)
            endswith(file, ".json") || continue
            d_json = open(JSON.parse, file)
            d_jl = repr(d_json)
            write(splitext(file)[1] * ".jl", d_jl)
        end
    end
end


text:
fn = "examples.md"
isfile(fn) && rm(fn)

for i = 1:6
    h = """
        Example $i
        ---------

        ![Example $i](http://www.cita.utoronto.ca/~nolta/julia/winston/examples/example$i.png)

        ``` julia"""
    run(`echo $h` >> fn)
    run(`cat ../examples/example$i.jl` >> fn)
    f = """
        ```
        """
    run(`echo $f` >> fn)
end


text:
@doc raw"""
    meanSquaredError(M, p, q)

Compute the (mean) squared error between the two
points `p` and `q` on the (power) manifold `M`.
"""
function meanSquaredError(M::mT, p, q) where {mT<:Manifold}
    return distance(M, p, q)^2
end
function meanSquaredError(M::PowerManifold, x, y)
    return 1 / prod(power_dimensions(M)) * sum(distance.(Ref(M.manifold), x, y) .^ 2)
end
@doc raw"""
    meanSquaredError(M,x,y)

Compute the (mean) squared error between the two
points `x` and `y` on the `PowerManifold` manifold `M`.
"""
function meanAverageError(M::Manifold, x, y)
    return distance(M, x, y)
end
function meanAverageError(M::PowerManifold, x, y)
    return 1 / prod(power_dimensions(M)) * sum(distance.(Ref(M.manifold), x, y))
end


text:
"""
    gmesh2fem(filename)

Extract points and elements from Gmesh file.
"""
function gmesh2fem(filename)
    error("Not implemented")
end


text:
using ParameterizedFunctions, PyDSTool, Test

f = @ode_def_bare Calcium begin
  dv = ( i + gl * (vl - v) - gca * 0.5 * (1 + tanh( (v-v1)/v2 )) * (v-vca) )/c
  dw = v-w
end vl vca i gl gca c v1 v2
p = [-60,120,0.0,2,4,20,-1.2,18]
u0 = [0;0]
tspan = [0;30]

dsargs = build_ode(f,u0,tspan,p)
#Solve the ODE
d = solve_ode(dsargs)
#using Plots
#plot(d[:t],d[:v])

ode = ds[:Generator][:Vode_ODEsystem](dsargs)
ode[:set](pars = Dict("i"=>-220))
ode[:set](ics  = Dict("v"=>-170))
PC = ds[:ContClass](ode)

bif = bifurcation_curve(PC,"EP-C",["i"],
                        max_num_points=450,
                        max_stepsize=2,min_stepsize=1e-5,
                        stepsize=2e-2,loc_bif_points="all",
                        save_eigen=true,name="EQ1",
                        print_info=true,calc_stab=true)

@test length(bif.changes) == 2

#=
using Plots
plot(bif,(:i,:v))
=#


text:
"""
    getCalmarRatio_percentage
calculates the calmarRatio for a given "profit and loss- series
Calmar Ratio = Average Annual Rate of Return / Maximum Drawdown
"""
function getCalmarRatio_percentage(profitAndLossArray::Array{Float64,1}, capitalAtBeginn::Float64, beginDate::Date, endDate::Date)::Float64
    averageAnnualReturnRate::Float64 = getAnnulizedReturnRateBetweenDatesYearly(profitAndLossArray,capitalAtBeginn,beginDate,endDate)
    maxDD_absolute::Float64 = getMaxDD(profitAndLossArray)
    maxDD_percentage = getMaxDDPercentage(maxDD_absolute,capitalAtBeginn)
    averageAnnualReturnRate_percentage = convertRateToPercent(averageAnnualReturnRate)
    return averageAnnualReturnRate_percentage/maxDD_percentage
end

max_stars_repo_name: RainerNeu/tradeStatistics

max_issues_repo_name: RainerNeu/tradeStatistics

max_forks_repo_name: RainerNeu/tradeStatistics


text:
## hack attempts to provide some of SymPy's simplify routines.
## could be seriously improved

## + expand  (built in)
##
## + expand_log
## + expand_power_exp
## + expand_power_base
## + expand_trig

## + logcombine (logsimp)
## + powsimp
## + powdenest
## * combsimp
## + trigsimp

## * factor (over Q)
## * collect
## * coeffs
## * cancel
## * apart

### --------------------------------------------------

export trig_simp, expand_trig

get_fun(ex::BasicType{Val{:Add}}) = :(+)
get_fun(ex::BasicType{Val{:Sub}}) = :(-)
get_fun(ex::BasicType{Val{:Mul}}) = :(*)
get_fun(ex::BasicType{Val{:Div}}) = :(/)
get_fun(ex::BasicType{Val{:Pow}}) = :(^)
#for t in (:Symbol, :Integer, :Rational, :Complex) ...
get_fun(ex::BasicType{Val{:Symbol}}) = :identity
get_fun(ex::BasicType{Val{:Integer}}) = :identity
get_fun(ex::BasicType{Val{:Rational}}) = :identity
get_fun(ex::BasicType{Val{:Complex}}) = :identity

function get_fun(u::BasicType)
    ex = Basic(u)
    as = args(ex)
    fn = head(ex) |> string |> lowercase |> Symbol
    fn
end

abstract type AbstractSimplify end

## identities

## rewrite product of sin and cos

## sin(x)*cos(y) = 1/2 * (sin(x+y) + sin(x-y)) ...
function sc(ex)
    @vars u_ v_ x___
    out = replace(ex,
                  sin(u_)*cos(v_)*x___,
                  (isempty(x___)?Basic(1):prod(x___)) * (sin(expand(v_ + u_)) + sin(expand(u_ - v_))))
    expand(out)
end

function ss(ex)
    @vars u_ v_ x___
    out = replace(ex,
                  sin(u_)*sin(v_)*x___,
                  x___*(1//2)*(cos(expand(u_ - v_)) - cos(expand(u_ + v_))))
    expand(out)
end

function cc(ex)
    @vars u_ v_ x___
    out = replace(ex,
                  cos(u_)*cos(v_)*x___,
                  x___ * (1//2)*(cos(expand(u_ + v_)) + cos(expand(u_ - v_))))
    expand(out)
end

## hyperbolic
function shch(ex)
    @vars u_ v_ x___
    out = replace(ex,
                  sinh(u_)*cosh(v_)*x___,
                  x___*(1//2)*(sinh(expand(v_ + u_)) + sinh(expand(u_ - v_))))
    expand(out)
end

function shsh(ex)
    @vars u_ v_ x___
    out = replace(ex,
                  sinh(u_)*sinh(v_)*x___,
                  x___*(1//2)*(cosh(expand(u_ + v_))- cosh(expand(u_ - v_)) ))
    expand(out)
end

function chch(ex)
    @vars u_ v_ x___
    out = replace(ex,
                  cosh(u_)*cosh(v_)*x___,
                  x___ * (1//2)*(cosh(expand(u_ + v_)) + cosh(expand(u_ - v_))))
    expand(out)
end

### expand trig -- sin(x+y) -> sin(x)*sin(y) + cos(x)*cos(y)
function sxy(ex)
    @vars x_ y_ x___
    ex = replace(ex, sin(x_ + y_) * x___ => x___ * (sin(x_)*cos(y_) + sin(y_) * cos(x_)))
#    m = match(ex, sin(x_*y_))
#    !ismatch(m) && return ex
    #    is XX Number(x_) or Number(y_) and integer and agreater 2 then write as 2x + (n-2)x and expand#...
    ex
end

function cxy(ex)
    @vars x_ y_ x___
    replace(ex, cos(x_ + y_) * x___ => x___ * (cos(x_)*cos(y_) - sin(x_) * sin(y_)))
end

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

##
## The basic _simp functions walk the syntax tree
## the specialized types then define what happens
## when the walking occurs with overrides as needed
function _simp(T::AbstractSimplify, ex::BasicType{Val{:Add}})
    sum(_simp(T, a) for a in args(Basic(ex)))
end

function _simp(T::AbstractSimplify, ex::BasicType{Val{:Mul}})
    prod(_simp(T, a) for a in args(Basic(ex)))
end

## "canonicalize" trig, hyperbolic, to use sin,cos or sinh, cosh
## Could make this a separate simplification, but that would
## involve walking the syntax tree again
_simp(::AbstractSimplify, ex::BasicType{Val{:Tan}}) = (u = args(Basic(ex));sin(u...)/cos(u...))
_simp(::AbstractSimplify, ex::BasicType{Val{:Cot}}) = (u = args(Basic(ex));cos(u...)/sin(u...))
_simp(::AbstractSimplify, ex::BasicType{Val{:Sec}}) = (u = args(Basic(ex));inv(cos(u...)))
_simp(::AbstractSimplify, ex::BasicType{Val{:Csc}}) = (u = args(Basic(ex));inv(sin(u...)))

_simp(::AbstractSimplify, ex::BasicType{Val{:Tanh}}) = (u = args(Basic(ex));sinh(u...)/cosh(u...))
_simp(::AbstractSimplify, ex::BasicType{Val{:Coth}}) = (u = args(Basic(ex));cosh(u...)/sinh(u...))
_simp(::AbstractSimplify, ex::BasicType{Val{:Sech}}) = (u = args(Basic(ex));inv(cosh(u...)))
_simp(::AbstractSimplify, ex::BasicType{Val{:Csch}}) = (u = args(Basic(ex));inv(sinh(u...)))

## leaves -- when args(ex) = []
function _simp(::AbstractSimplify, ex::T) where {T <: SymEngine.BasicNumber}
    Basic(ex)
end

function _simp(::AbstractSimplify, ex::BasicType{Val{:Symbol}})
    Basic(ex)
end

## catch all
function _simp(T::AbstractSimplify, ex::BasicType)
#    println("_simp: $ex $(typeof(ex))")
    eval(Expr(:call, get_fun(ex), broadcast(_simp, T, args(Basic(ex)))...))
end
_simp(T::AbstractSimplify, ex::Basic) = _simp(T, BasicType(ex))

## iterate until a fixed point
function _fixed_point(T::AbstractSimplify, ex, n=20)
    for i in 1:n
        ex1 = _simp(T, BasicType(ex))
        ex1 == ex && return expand(ex)
        ex = ex1
    end
    expand(ex)
end

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

## check that walking works
struct BlankSimp <: AbstractSimplify end

##################################################
## TrigSimp
struct TrigSimp <: AbstractSimplify end

"""
    trig_simp

Expand trig products such as `sin(x)*cos(y)` into sums using double angle
formulas
"""
trig_simp(ex::Basic) = _fixed_point(TrigSimp(), ex)

function _simp(T::TrigSimp, ex::BasicType{Val{:Pow}})
    a,b = args(Basic(ex))
    if  head(b) == :Integer && b > 1
        @vars u_
        op = head(a)
        op == :Sin  && return replace(a^2, sin(u_)^2,   1//2*(1 - cos(2*u_)))  * a^(b-2)
        op == :Cos  && return replace(a^2, cos(u_)^2,   1//2*(1 + cos(2*u_)))  * a^(b-2)
        op == :Sinh && return replace(a^2, sinh(u_)^2, -1//2*(1 - cosh(2*u_))) * a^(b-2)
        op == :Cosh && return replace(a^2, cosh(u_)^2,  1//2*(1 + cosh(2*u_))) * a^(b-2)
    elseif head(b) == :Integer && b <= -1
        return _simp(T, a)^b
    end
    return ex
end

## products apply the simplification rules for sin*cos,...
function _simp(T::TrigSimp, ex::BasicType{Val{:Mul}})
    v = Basic(ex)
    u = v |> cc |> ss |> sc |> chch |> shsh |> shch
    u = expand(u)
    head(u) != :Mul && return u
    out = Basic(1)

    ## out =prod(_simp(T, a) for a in args(Basic(u)))
    for a in args(Basic(u))
        out = out * _simp(T, a)
    end
    return out
end

##################################################
## expand trig
struct ExpandTrig <: AbstractSimplify end

"""
    expand_trig

Apply double angle formulas to express as products of sines and cosines
"""
expand_trig(ex::Basic) = _fixed_point(ExpandTrig(), ex)

_simp(T::ExpandTrig, ex::BasicType{Val{:Sin}}) = sxy(Basic(ex))
_simp(T::ExpandTrig, ex::BasicType{Val{:Cos}}) = cxy(Basic(ex))

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


text:
# by lazarusA # HIDE
#using GLMakie, Random # HIDE
using CairoMakie, Random
CairoMakie.activate!() #HIDE
let
    Random.seed!(123)
    n = 15
    x, y, color = rand(n), rand(n), rand(n)
    cmaps=[:cool,:viridis,:plasma,:inferno,:thermal,:leonardo,:winter,:spring,:ice]

    function FigGridHeatSharedCbar()
        fig = Figure(resolution = (700, 500))
        c = 1
        for i in 1:2, j in 1:2:3
            ax=Axis(fig[i, j],aspect = 1,xgridvisible = false,ygridvisible = false)
            pnts = heatmap!(rand(10,10), colormap= cmaps[j], colorrange=(0, 1))
            ax.xticks = [1,10]
            ax.yticks = [1,10]
            ax.xticklabelsize = 16
            ax.yticklabelsize = 16

            cbar = Colorbar(fig, pnts, ticklabelsize = 14, height = Relative(3/4),
            tickwidth = 2)
            cbar.ticks = [0, 0.5, 1]
            fig[1:2, j+1] = cbar
            c+=1
        end
        fig
    end
    fig = FigGridHeatSharedCbar()
    fig
    save(joinpath(@__DIR__, "output", "heatmap2SharedCbars.png"), fig, px_per_unit = 2.0) # HIDE
end
using Pkg # HIDE
Pkg.status(["CairoMakie", "Random"]) # HIDE

max_stars_repo_name: MisterBiggs/BeautifulMakie

max_issues_repo_name: MisterBiggs/BeautifulMakie

max_forks_repo_name: MisterBiggs/BeautifulMakie


text:
using BenchmarkTools

# generic tools for evaluation
"""
    testfuns(F, args...; kwargs...)

Benchmark the functions in F using `args` and `kwargs`.
"""
function testfuns(F::Union{Array{Function,1}, Array{Type,1}}, args...; kwargs...)
    # set up benchmark suite
    suite = BenchmarkGroup()
    for f in F
        suite[string(f)] = @benchmarkable $(f)($(args)...; $(kwargs)...)
    end
    # tune and run suite
    tune!(suite)
    results = run(suite, verbose=get(kwargs, :verbose, false))
    # show results
    for f in F
        println("\nResult for $(f):\n$("-"^length("Result for $(f):"))")
        display(results[string(f)])
    end
    medians = sort!([(fun, median(trial.times)) for (fun, trial) in results], by=x->x[2])
    speedup = round(Int,medians[end][2] / medians[1][2])
    println("\n$(speedup)x speedup (fastest: $(medians[1][1]), slowest: $(medians[end][1]))\n")
    return
end

"""
    equality(F::Array{Function,1}, args...; kwargs...)

Verify the result of an array of functions F is the same using `args` and `kwargs`.

### Notes
The first function in the array is used as the reference
"""
function equality(F::Array{Function,1}, args...; kwargs...)
    return all([equality(F[1], F[i], args...; kwargs...) for i in 2:length(F)])
end

"""
    equality(f::Function, g::Function, args...; kwargs...)

Verify the result of function f and g is the same (element-wise) using `args` and `kwargs`.
"""
function equality(f::Function, g::Function, args...; kwargs...)
    return all(isequal.(map(f->f(args...; kwargs...), [f, g])...))
end

# avoid globals (general rule unless constants, in that case define them as such)

"""
    renew(char="", n::Int=0)

Slow implementation of string joining
"""
function renew(char="", n::Int=0; kwargs...)
    res = ""
    for i in 1:n
        res *= char
    end

    return res
end

"""
    buffer(char="", n::Int=0)

Fast implementation of string joining
"""
function buffer(char="", n::Int=0; kwargs...)
    res = IOBuffer()
    for i in 1:n
        print(res, char)
    end

    return String(take!(res))
end

"""
    demo1

- Illustration of implementation issues

speedup ∝ length
"""
function demo1()
    F = [renew, buffer];
    args=("a", 10000)
    # test validity
    @info "Testing demo 1"
    @assert equality(F, args...)
    # Run benchmark
    @info "Benching demo 1"
    testfuns(F, args...)
end

function f_alloc(x::Vector{T}, y::Array{T,2}; props::Vector{T}, kwargs...) where T
    # read properties
    prop_A = props[1]
    prop_B = props[2]
    prop_C = props[3]
    prop_D = props[4]
    # initialise
    result = zeros(T, size(x))
    result .+= x
    # update result
    result = result .+ prop_A .* y[:,1]
    result = result .+ prop_B .* y[:,2]
    result = result .+ prop_C .* y[:,3]
    result = result .+ prop_D .* y[:,4]
    return result
end

struct Props{T<:AbstractFloat}
    prop_A::T
    prop_B::T
    prop_C::T
    prop_D::T
end

@views function f_noalloc(x::Vector{T}, y::Array{T,2}; p::Props, kwargs...) where T
    return x .+ p.prop_A .* y[:,1] .+ p.prop_B .* y[:,2] + p.prop_C .* y[:,3] + p.prop_D .* y[:,4]
end

"""
    demo2()

- Avoiding unnecessary allocations
- Using views
- Avoiding fields with abstract types
"""
function demo2()
    F = [f_alloc, f_noalloc]
    N = 10000
    x = rand(N)
    y = rand(N,4)
    args = [x, y]
    props = Float64[1,2,3,4]
    p = Props(props...)
    kwargs = Dict(:props => props, :p => p)
    # test validity
    @info "Testing demo 2"
    @assert equality(F, args...; kwargs...)
    # run benchmark
    @info "Benching demo 2"
    testfuns(F, args...; kwargs...)
end

# avoid changing type of a variable
function f_typechange(x)
    res = 0
    for val in x
        res += val
    end
    return res
end

function f_notypechange(x::Array{T,1}) where T <: Number
    res = zero(T)
    for val in x
        res += val
    end
    return res
end

"""
    demo3()

- avoid changing type of a variable
"""
function demo3()
    F = [f_typechange, f_notypechange]
    N = 2
    args = [rand(N)]
    # test validity
    @info "Testing demo 3"
    @assert equality(F, args...)
    # run benchmark
    @info "Benching demo 3"
    testfuns(F, args...)
end

f_norm(x) = 3x.^2 + 4x + 7x.^3;
f_vec(x) = @. 3x.^2 + 4x + 7x.^3;

"""
    demo 4

- use vectorisation where possible
"""
function demo4()
    F = [f_norm, f_vec]
    N::Int = 1e6
    args = [rand(N)]
    # test validity
    @info "Testing demo 4"
    @assert equality(F, args...)
    # run benchmark
    @info "Benching demo 4"
    testfuns(F, args...)
end

function f_nopre(x::Vector; kwargs...)
    y = similar(x)
    for i in eachindex(x)
        y[i] = sum(x[1:i])
    end
    return y
end

function f_pre!(x::Vector{T}; y::Vector{T}, kwargs...) where {T}
    for i in eachindex(x)
        y[i] = sum(x[1:i])
    end
    return y
end

function f_preboost!(x::Vector{T}; y::Vector{T}, kwargs...) where {T}
    Threads.@threads for i in eachindex(x)
        @inbounds y[i] = sum(x[1:i])
    end
    return y
end

@views function f_preturboboost!(x::Vector{T}; y::Vector{T}, kwargs...) where {T}
    Threads.@threads for i in eachindex(x)
        @inbounds y[i] = sum(x[1:i])
    end
    return y
end

"""
    demo5()

- preallocate outputs if possible
- make use of threads if possible
- make use of views if possible
- skip index check when sure of dimensions
"""
function demo5()
    F = [f_nopre, f_pre!, f_preboost!, f_preturboboost!]
    @info "Currently running on $(Threads.nthreads()) threads"
    N::Int = 1e4
    x::Vector = rand(N)
    y_pre = similar(x)
    args = [x]
    kwargs = Dict( :y => y_pre)
    # test validity
    @info "Testing demo 5"
    @assert equality(F, args...; kwargs...)
    # run benchmark
    @info "Benching demo 5"
    testfuns(F, args...; kwargs...)
end

function grow(N::Int; kwargs...)
    res = []
    for _ in 1:N
        push!(res, 1)
    end
    return res
end

function nogrow(N::Int; x::Array{T,1}) where{T}
    res = similar(x)
    for i in 1:N
        res[i] = 1
    end
    return res
end

"""
    demo6()

- growing vs. allocating for known dimensions
"""
function demo6()
    N = 10000
    x = rand(N)
    F = [grow, nogrow]
    args = [N]
    kwargs = Dict(:x => x)
    # test validity
    @info "Testing demo 6"
    @assert equality(F, args...; kwargs...)
    # run benchmark
    @info "Benching demo 6"
    testfuns(F, args...; kwargs...)
end

struct mystruct
    a
    b
    c
end

struct mytypedstruct{T}
    a::T
    b::T
    c::T
end

mutable struct mymutablestruct
    a
    b
    c
end

"""
    demo7

- creation time of mutable vs unmutable structs and the parameter version
"""
function demo7()
    structs = [mystruct, mytypedstruct, mymutablestruct]
    @info "Benching demo 7"
    testfuns(structs, [1;2;3]...)
end

function reader(s::mystruct)
    return (s.a, s.b, s.c)
end

function reader(s::mytypedstruct)
    return (s.a, s.b, s.c)
end

function reader(sm::mymutablestruct)
    return (sm.a, sm.b, sm.c)
end

"""
    demo8

- read times of fields of mutable vs unmutable structs and the parameter version
"""
function demo8()
    vals = [1;2;3]
    @info "Benching demo8"
    for arg in [mystruct, mytypedstruct, mymutablestruct]
        println("\nResult for $(arg):\n$("-"^length("Result for $(arg):"))")
        x = arg(vals...)
        display(@benchmark reader($(x)))
    end
end

function main(N::Int=8)
    for i in 1:N
        eval(Symbol("demo$(i)"))()
    end
end

main()


text:
@testset "leftgroupjoin" begin
    l = [1,2,3,4]
    r = [0,1,2]
    ans1 = dictionary([1 => [(1,1)],
                       2 => [(2,2)],
                       3 => [],
                       4 => []])
    ans2 = dictionary([false => [(1, 1), (3, 1)],
                       true => [(2, 0), (2, 2), (4, 0), (4, 2)]])

    @test leftgroupjoin(identity, identity, tuple, l ,r) == ans1
    @test leftgroupjoin(iseven, iseven, tuple, l,  r) == ans2
    @test leftgroupjoin(iseven, iseven, tuple, ==, l, r) == ans2
end


text:
using ROCArrays
using HSARuntime
using LinearAlgebra
using Test

include("util.jl")

# copy-pasta from GPUArrays/src/testsuite.jl
convert_array(f, x) = f(x)
convert_array(f, x::Base.RefValue) = x[]
function compare(f, AT::Type{ROCArray}, xs...; kwargs...)
    cpu_in = convert_array.(copy, xs)
    gpu_in = convert_array.(AT, xs)
    cpu_out = f(cpu_in...; kwargs...)
    gpu_out = f(gpu_in...; kwargs...)
    collect(cpu_out) ≈ collect(gpu_out)
end
# copy-pasta from CuArrays/test/runtests.jl
testf(f, xs...; kwargs...) = compare(f, ROCArray, xs...; kwargs...)

const agent = HSARuntime.get_default_agent()
agent_name = HSARuntime.get_name(agent)
agent_isa = HSARuntime.get_first_isa(agent)
@info "Testing using device $agent_name with ISA $agent_isa"

@testset "ROCArrays" begin

@testset "ROCArrays External Libraries" begin
    isdefined(ROCArrays, :rocBLAS) ? include("blas.jl") : @test_skip "rocBLAS"
    isdefined(ROCArrays, :rocFFT) ? include("fft.jl") : @test_skip "rocFFT"
end

end # @testset "ROCArrays"


text:
# Copyright (c) 2018: Matthew Wilhelm & Matthew Stuber.
# This code is licensed under MIT license (see LICENSE.md for full details)
#############################################################################
# EAGO
# A development environment for robust and global optimization
# See https://github.com/PSORLab/EAGO.jl
#############################################################################
# Defines the EmptyNLPEvaluator needed to provide a constructor for
# the EAGO.Optimizer prior to loading JuMP's NLPEvaluator.
#############################################################################

struct EmptyNLPEvaluator <: MOI.AbstractNLPEvaluator
    _current_node::NodeBB
    has_nlobj::Bool
end
EmptyNLPEvaluator() = EmptyNLPEvaluator(NodeBB(),false)
set_current_node!(x::EmptyNLPEvaluator, n::NodeBB) = ()

MOI.features_available(::EmptyNLPEvaluator) = [:Grad, :Jac, :Hess]
MOI.initialize(::EmptyNLPEvaluator, features) = nothing
MOI.eval_objective(::EmptyNLPEvaluator, x) = NaN
function MOI.eval_constraint(::EmptyNLPEvaluator, g, x)
    @assert length(g) == 0
    return
end
MOI.eval_objective_gradient(::EmptyNLPEvaluator, g, x) = nothing
MOI.jacobian_structure(::EmptyNLPEvaluator) = Tuple{Int64,Int64}[]
MOI.hessian_lagrangian_structure(::EmptyNLPEvaluator) = Tuple{Int64,Int64}[]
function MOI.eval_constraint_jacobian(::EmptyNLPEvaluator, J, x)
    @assert length(J) == 0
    return
end
function MOI.eval_hessian_lagrangian(::EmptyNLPEvaluator, H, x, σ, μ)
    @assert length(H) == 0
    return
end

empty_nlp_data() = MOI.NLPBlockData([], EmptyNLPEvaluator(), false)


text:
# functions to update incycle, containRoot, gammaz
# Claudia March 2015
#####################

# --------------------------------------- update incycle, root, gammaz -------------------------------------------

# function to update inCycle (with priority queue) after becoming part of a network
# based on program 3 CS367 with priority queue
# expected to be much faster than the other two udpateInCycle (queue and recursive)
# input: hybrid node around which we want to update inCycle
# returns tuple: flag, nocycle, array of edges changed, array of nodes changed
#   flag: false if cycle intersects existing cycle or number of nodes in cycle < 3
#         (there is the possibility of returning edges in intersection: path)
#         true if cycle does not intersect existing cycle
#   nocycle: true if there is no cycle (instead of error). it is used in addHybridization
# calculates also the number of nodes in the cycle and put as hybrid node attribute "k"
# warning: it is not checking if hybrid node or minor hybrid edge
#          were already part of a cycle (inCycle!= -1)
#          But it is checking so for the other edges in cycle
# warning: it needs extra things: visited attribute, prev attribute
#          unlike updateInCycle recursive, but it is expected
#          to be much faster
function updateInCycle!(net::HybridNetwork,node::Node)
    if(node.hybrid)
        start = node;
        node.inCycle = node.number;
        node.k = 1;
        hybedge = getHybridEdge(node);
        hybedge.inCycle = node.number;
        last = getOtherNode(hybedge,node);
        dist = 0;
        queue = PriorityQueue();
        path = Node[];
        net.edges_changed = Edge[];
        net.nodes_changed = Node[];
        push!(net.edges_changed,hybedge);
        push!(net.nodes_changed,node);
        found = false;
        net.visited = [false for i = 1:size(net.node,1)];
        enqueue!(queue,node,dist);
#        println("start is $(start.number), last is $(last.number)")
        while(!found)
            if(isempty(queue))
                return false, true, net.edges_changed, net.nodes_changed
            else
                #println("at this moment, queue is $([n.number for n in queue])")
                curr = dequeue!(queue);
 #               println("curr $(curr.number)")
                if(isEqual(curr,last))
 #                   println("encuentra a la last $(curr.number)")
                    found = true;
                    push!(path,curr);
                else
                    if(!net.visited[getIndex(curr,net)])
 #                       println("curr not visited $(curr.number)")
                        net.visited[getIndex(curr,net)] = true;
                        if(isEqual(curr,start))
                            #println("curr is start")
                            for e in curr.edge
                                if(!e.hybrid || e.isMajor)
                                    other = getOtherNode(e,curr);
 #                                   println("other is $(other.number), pushed to queue")
                                    other.prev = curr;
                                    dist = dist+1;
                                    enqueue!(queue,other,dist);
                                end
                            end
                        else
                            for e in curr.edge
                                if(!e.hybrid || e.isMajor)
                                    other = getOtherNode(e,curr);
 #                                   println("other is $(other.number)")
                                    if(!other.leaf && !net.visited[getIndex(other,net)])
 #                                       println("dice que other $(other.number) no es leaf ni visited, lo mete a queue")
                                        other.prev = curr;
                                        dist = dist+1;
                                        enqueue!(queue,other,dist);
                                    end
                                end
                            end
                        end
                    end
                end
            end
        end # end while
 #       println("after while, path has $([n.number for n in path])")
        curr = pop!(path);
 #       println("path should be empty here: $(isempty(path))")
        while(!isEqual(curr, start))
 #           println("curr es $(curr.number)")
            if(curr.inCycle!= -1)
                push!(path,curr);
                curr = curr.prev;
            else
                curr.inCycle = start.number;
                push!(net.nodes_changed, curr);
                node.k  =  node.k + 1;
                edge = getConnectingEdge(curr,curr.prev);
                edge.inCycle = start.number;
                push!(net.edges_changed, edge);
                curr = curr.prev;
            end
        end
        if(!isempty(path) || node.k<3)
            @debug "warning: new cycle intersects existing cycle"
            return false, false, net.edges_changed, net.nodes_changed
        else
            return true, false, net.edges_changed, net.nodes_changed
        end
    else
        error("node is not hybrid")
    end
end

"""
    updateContainRoot!(HybridNetwork, Node)
    traverseContainRoot!(Node, Edge, edges_changed::Array{Edge,1}, rightDir::Vector{Bool})

The input `node` to `updateContainRoot!` must be a hybrid node
(can come from searchHybridNode).
`updateContainRoot!` starts at the input node and calls `traverseContainRoot!`,
which traverses the network recursively.
By default, containRoot attributes of edges are true.
Changes `containRoot` to false for all the visited edges: those
below the input node, but not beyond any other hybrid node.

`updateContainRoot!` Returns a `flag` and an array of edges whose
containRoot has been changed from true to false.
`flag` is false if the set of edges to place the root is empty

In `traverseContainRoot!`, `rightDir` turns false if hybridizations
have incompatible directions (vector of length 1, to be modified).

Warning:

- does *not* update `containRoot` of minor hybrid edges.
- assumes correct `isMajor` attributes: to stop the recursion at minor hybrid edges.
- assumes correct hybrid attributes of both nodes & edges: to check if various
  hybridizations have compatible directions.
  For each hybrid node that is encountered, checks if it was reached
  via a hybrid edge (ok) or tree edge (not ok).

`rightDir`: vector of length 1 boolean, to be mutable and modified by the function
"""
function traverseContainRoot!(node::Node, edge::Edge, edges_changed::Array{Edge,1}, rightDir::Vector{Bool})
    if node.hybrid
        if edge.hybrid
            edge.isMajor || error("hybrid edge $(edge.number) is minor and we should not traverse the graph through minor edges")
            DEBUGC && @debug "traverseContainRoot reaches hybrid node $(node.number) through major hybrid edge $(edge.number)"
            rightDir[1] &= true  # This line has no effect: x && true = x
        else #approach hybrid node through tree edge => wrong direction
            rightDir[1] &= false # same as rightDir[1] = false: x && false = false
            DEBUGC && @debug "traverseContainRoot reaches hybrid node $(node.number) through tree edge $(edge.number), so rightDir $(rightDir[1])"
        end
    elseif !node.leaf
        for e in node.edge
            if !isEqual(edge,e) && e.isMajor # minor edges avoided-> their containRoot not updated
                other = getOtherNode(e,node);
                if e.containRoot # only considered changed those that were true and not hybrid
                    DEBUGC && @debug "traverseContainRoot changing edge $(e.number) to false, at this moment, rightDir is $(rightDir[1])"
                    e.containRoot = false;
                    push!(edges_changed, e);
                end
                traverseContainRoot!(other,e, edges_changed, rightDir);
            end
        end
    end
end

# node: must be hybrid node (can come from searchHybridNode)
# return flag, array of edges changed
#        flag: false if the set of edges to place the root is empty
@doc (@doc traverseContainRoot!) updateContainRoot!
function updateContainRoot!(net::HybridNetwork, node::Node)
    node.hybrid || error("node $(node.number )is not hybrid, cannot update containRoot")
    net.edges_changed = Edge[];
    rightDir = [true] #assume good direction, only changed if found hybrid node through tree edge
    for e in node.edge
        if !e.hybrid
            other = getOtherNode(e,node);
            e.containRoot = false;
            push!(net.edges_changed,e);
            traverseContainRoot!(other,e, net.edges_changed,rightDir);
        end
    end
    if !rightDir[1] || all((e->!e.containRoot), net.edge)
        return false,net.edges_changed
    else
        return true,net.edges_changed
    end
end

# function to identify if the network is one of the pathological cases
# see ipad notes: k = 0 (nonidentifiable), k = 1 (nonidentifiable
# "extreme bad triangle", "very" bad triangle I,II) k = 2 (bad diamond
# I,II) also checks if hybrid node has leaf child, in which case,
# major edge is non identifiable
# input: hybrid node around which to check (can come from searchHybridNode)
#        updates gammaz with whatever
# edge lengths are originally in the network
#        allow = true, returns true always, used when reading topology
# returns net.hasVeryBadTriangle, array of edges changed (istIdentifiable, except hybrid edges)
#         false if the network has extremely/very bad triangles
# warning: needs to have updateInCycle already done as it needs inCycle, and k
# check: assume any tree node that has hybrid Edge has only
# one tree edge in cycle (true?)
# updates net.numBad attribute when found a bad diamond I
function updateGammaz!(net::HybridNetwork, node::Node, allow::Bool)
    node.hybrid || error("node $(node.number) is not hybrid, cannot updategammaz")
    node.k != -1 || error("update in cycle should have been run before: node.k not -1")
    node.isExtBadTriangle = false
    node.isVeryBadTriangle = false
    node.isBadTriangle = false
    node.isBadDiamondI = false
    node.isBadDiamondII = false
    net.edges_changed = Edge[];
    edge_maj, edge_min, tree_edge2 = hybridEdges(node);
    other_maj = getOtherNode(edge_maj,node);
    other_min = getOtherNode(edge_min,node);
    node.k > 2 || return false, []
    if(node.k == 4) # could be bad diamond I,II
#        net.numTaxa >= 5 || return false, [] #checked inside optTopRuns now
        edgebla,edge_min2,tree_edge3 = hybridEdges(other_min);
        edgebla,edge_maj2,tree_edge1 = hybridEdges(other_maj);
        other_min2 = getOtherNode(edge_min2,other_min);
        isLeaf1 = getOtherNode(tree_edge1,other_maj);
        isLeaf2 = getOtherNode(tree_edge2,node);
        isLeaf3 = getOtherNode(tree_edge3,other_min);
        tree_edge4 = nothing;
        for e in other_min2.edge
            if(isa(tree_edge4,Nothing) && e.inCycle == -1 && !e.hybrid)
                tree_edge4 = e;
            end
        end
        if(isEqual(other_min2,getOtherNode(edge_maj2,other_maj)) && isLeaf1.leaf && isLeaf2.leaf && isLeaf3.leaf) # bad diamond I
            @debug "bad diamond I found"
            net.numBad += 1
            node.isBadDiamondI = true;
            other_min.gammaz = edge_min.gamma*edge_min2.z;
            other_maj.gammaz = edge_maj.gamma*edge_maj2.z;
            edge_min2.istIdentifiable = false;
            edge_maj2.istIdentifiable = false;
            edge_maj.istIdentifiable = false;
            edge_min.istIdentifiable = false;
            push!(net.edges_changed,edge_min2);
            push!(net.edges_changed,edge_min);
            push!(net.edges_changed,edge_maj2);
            push!(net.edges_changed,edge_maj);
        elseif(isEqual(other_min2,getOtherNode(edge_maj2,other_maj)) && isLeaf1.leaf && !isLeaf2.leaf && isLeaf3.leaf && getOtherNode(tree_edge4,other_min2).leaf) # bad diamond II
            @debug "bad diamond II found"
            node.isBadDiamondII = true;
            setLength!(edge_maj,edge_maj.length+tree_edge2.length)
            setLength!(tree_edge2,0.0)
            push!(net.edges_changed,tree_edge2)
            tree_edge2.istIdentifiable = false
            edge_maj.istIdentifiable = true
            edge_min.istIdentifiable = true
        end
    elseif(node.k == 3) # could be extreme/very bad triangle or just bad triangle
        if(net.numTaxa <= 5)
            @debug "extremely or very bad triangle found"
            node.isVeryBadTriangle = true
            net.hasVeryBadTriangle = true
        elseif(net.numTaxa >= 6)
            edgebla,tree_edge_incycle,tree_edge1 = hybridEdges(other_min);
            edgebla,edgebla,tree_edge3 = hybridEdges(other_maj);
            isLeaf1 = getOtherNode(tree_edge1,other_min);
            isLeaf2 = getOtherNode(tree_edge2,node);
            isLeaf3 = getOtherNode(tree_edge3,other_maj);
            if isLeaf1.leaf || isLeaf2.leaf || isLeaf3.leaf
                nl = count([l.leaf for l in [isLeaf1,isLeaf2,isLeaf3]])
                if nl >= 2
                    @debug "warning: extremely bad triangle found"
                    node.isExtBadTriangle = true;
                    net.hasVeryBadTriangle = true
                elseif nl == 1
                    @debug "warning: bad triangle I or II found"
                    node.isVeryBadTriangle = true;
                    net.hasVeryBadTriangle = true
                end
            else
                node.isBadTriangle = true
                setLength!(edge_maj,edge_maj.length+tree_edge2.length)
                setLength!(tree_edge2,0.0)
                tree_edge2.istIdentifiable = false
                push!(net.edges_changed, tree_edge2);
            end
        end
    end #ends the search for bad things
    if(node.k > 3 && !node.isBadDiamondI && !node.isBadDiamondII)
        #println("si entra el ultimo if de k>3 y no bad diamondI,II")
        edgebla,tree_edge_incycle,tree_edge1 = hybridEdges(other_min);
        if(!tree_edge_incycle.istIdentifiable)
            tree_edge_incycle.istIdentifiable = true;
            push!(net.edges_changed,tree_edge_incycle);
        end
        edge_maj.istIdentifiable = isEdgeIdentifiable(edge_maj)
        edge_min.istIdentifiable = isEdgeIdentifiable(edge_min)
    end
    isBadTriangle(node) == net.hasVeryBadTriangle || error("node $(node.number) is very bad triangle but net.hasVeryBadTriangle is $(net.hasVeryBadTriangle)")
    if(allow)
        return true, net.edges_changed
    else
        return !net.hasVeryBadTriangle, net.edges_changed
    end
end

updateGammaz!(net::HybridNetwork, node::Node) = updateGammaz!(net, node, false)

#function to check if edge should be identifiable
#it is not only if followed by leaf, or if a newly converted tree edge
function isEdgeIdentifiable(edge::Edge)
    if(edge.hybrid)
        node = edge.node[edge.isChild1 ? 1 : 2]
        #println("is edge $(edge.number) identifiable, node $(node.number)")
        node.hybrid || error("hybrid edge $(edge.number) pointing at tree node $(node.number)")
        major,minor,tree = hybridEdges(node)
        #println("major $(major.number), minor $(minor.number), tree $(tree.number)")
        if(getOtherNode(tree,node).leaf)
            return false
        else
            return true
        end
    else
        if(reduce(&,[!edge.node[1].leaf,!edge.node[2].leaf]))
            if(!edge.node[1].hybrid && !edge.node[2].hybrid && !edge.fromBadDiamondI)
                return true
            elseif(edge.node[1].hybrid || edge.node[2].hybrid)
                ind = edge.node[1].hybrid ? 1 : 2
                if(!edge.node[ind].isBadDiamondII && !edge.node[ind].isBadTriangle)
                    return true
                else
                    return false
                end
            end
        else
            return false
        end
    end
end

# function to update the net.partition attribute along a cycle formed
# by nodesChanged vector (obtained from updateInCycle)
# warning: needs updateInCycle for all hybrids before running this
function updatePartition!(net::HybridNetwork, nodesChanged::Vector{Node})
    length(nodesChanged) > 2 || error("incycle with only 2 nodes in it after updateGammaz")
    #println("nodesChanged are $([n.number for n in nodesChanged])")
    if(net.numHybrids == 0)
        net.partition = Partition[]
    end
    for n in nodesChanged
        if(length(n.edge) == 3) #because we are allowing the root to have only two edges when read from parenthetical format
            edge = nothing
            for e in n.edge
                if(e.inCycle == -1)
                    edge = e
                end
            end
            !isa(edge,Nothing) || error("one edge in n.edge for node $(n.number) should not be in cycle")
            descendants = [edge]
            cycleNum = [nodesChanged[1].inCycle]
            getDescendants!(getOtherNode(edge,n),edge,descendants,cycleNum)
            !isempty(descendants) || error("descendants is empty for node $(n.number)")
            @debug "for node $(n.number), descendants are $([e.number for e in descendants]), and cycleNum is $(cycleNum)"
            partition = Partition(cycleNum,descendants)
            if(!isPartitionInNet(net,partition)) #need to check not already added by other hybrid nodes
                push!(net.partition, partition)
            end
        end
    end
end

function choosePartition(net::HybridNetwork)
    all((n->(length(n.edges) == 1)), net.partition) && return 0 #cannot put any hyb
    all((n->(length(n.edges) == 3)), net.partition) && return 0 #can only put very bad triangles
    partition = Int[] #good partitions
    for i in 1:length(net.partition)
        if(length(net.partition[i].edges) > 3)
            push!(partition,i)
        end
    end
    isempty(partition) && return 0
    length(partition) == 1 && return partition[1]
    index1 = round(Integer,rand()*size(partition,1));
    while(index1 == 0 || index1 > length(partition))
        index1 = round(Integer,rand()*size(partition,1));
    end
    @debug "chosen partition $([n.number for n in net.partition[partition[index1]].edges])"
    return partition[index1]
end

# based on getDescendants on readData.jl but with vector of edges, instead of nodes
# finds the partition corresponding to the node and edge in the cycle
# used in chooseEdgesGamma and to set net.partition
# cycleNum is a variable that will save another hybrid node number if found
function getDescendants!(node::Node, edge::Edge, descendants::Vector{Edge}, cycleNum::Vector{Int})
    @debug "getDescendants of node $(node.number) and edge $(edge.number)"
    if(node.inCycle != -1)
        push!(cycleNum,node.inCycle)
    elseif(!node.leaf && node.inCycle == -1)
        for e in node.edge
            if(!isEqual(edge,e) && e.isMajor)
                push!(descendants,e)
                getDescendants!(getOtherNode(e,node),e,descendants,cycleNum)
            end
        end
    end
end


text:
# Elementwise binary functions for arrays of the same size

import Base: .+, .-, .*, ./, .^, max, min, .==, .>, .>=, .<, .<=, +, -

cuda11 = [
("add",".+","xi+yi"),
("sub",".-","xi-yi"),
("mul",".*","xi*yi"),
("div","./","xi/yi"),
("pow",".^","pow(xi,yi)"),
("max","max","(xi>yi?xi:yi)"),
("min","min","(xi<yi?xi:yi)"),
("eq",".==","xi==yi"),
("gt",".>","xi>yi"),
("ge",".>=","xi>=yi"),
("lt",".<","xi<yi"),
("le",".<=","xi<=yi"),
# "hypot",
# "rhypot",
# "atan2",
# "frexp",
# "ldexp",
# "scalbn",
# "scalbln",
# "jn",
# "yn",
# "fmod",
# "remainder",
# "mod",
# "fdim",
("invxback","invxback","(-xi*yi*yi)"),
("reluback","reluback","(yi>0?xi:0)"),
("sigmback","sigmback","(xi*yi*(1-yi))"),
("tanhback","tanhback","(xi*(1-yi*yi))"),
]

# These are used in cuda12.jl as special cases of equal sized array,array ops.
# Here we'll just define some functions that specifically do not have broadcasting.

(+){T}(x::KnetArray{T},y::KnetArray{T})=(size(x)==size(y)||throw(DimensionMismatch("$(map(size,(x,y)))"));(.+)(x,y))
(-){T}(x::KnetArray{T},y::KnetArray{T})=(size(x)==size(y)||throw(DimensionMismatch("$(map(size,(x,y)))"));(.-)(x,y))
#(*){T}(x::KnetArray{T},y::KnetArray{T})=(.*)(x,y) # This is matmul
#(/){T}(x::KnetArray{T},y::KnetArray{T})=(./)(x,y) # This is another linalg op


text:
include("density_test.jl")
include("height_test.jl")
include("expansion_test.jl")


text:
function producer(c::Channel)
   put!(c, "start")
   for n=1:4
       put!(c, 2n)
   end
   put!(c, "stop")
end

ch = Channel(producer)
println(take!(ch))
println(take!(ch))
println(take!(ch))
println(take!(ch))
println(take!(ch))
println(take!(ch))
# error, can't access 'closed channel
# println(take!(ch))

max_stars_repo_name: mkanenobu/trashbox

max_issues_repo_name: mkanenobu/trashbox

max_forks_repo_name: mkanenobu/trashbox


text:
#Pkg.add("UnstructuredGrids")
using UnstructuredGrids

# Generate a grid from its nodal coordinates, cell connectivities,
# and cell types

# Define coordinates
coords = zeros(2,9)
coords[:,1] = [1,1]
coords[:,2] = [3,1]
coords[:,3] = [4,1]
coords[:,4] = [1,2]
coords[:,5] = [2,2]
coords[:,6] = [1,3]
coords[:,7] = [2,3]
coords[:,8] = [3,3]
coords[:,9] = [4,3]

# Define connectivity
connect = [1,2,5,4,2,3,9,4,5,7,6,5,8,7,5,2,8,2,9,8]
offsets = [1,      5,    8,      12,   15,   18,   21]

# Define cell types
using UnstructuredGrids.RefCellGallery: SQUARE, TRIANGLE
refcells = [SQUARE, TRIANGLE]
types = [1,2,1,2,2,2]

# Generate the UGrid object
grid = UGrid(connect,offsets,types,refcells,coords)

# Export grid into vtk format
#writevtk(grid,"grid")

# Generate a global numbering for the edges (1d objects) of the grid,
# and find which edges are on the boundary of each cell
n=1
cell_to_edges = generate_cell_to_faces(n,grid)

#@show cell_to_edges.list
#@show cell_to_edges.ptrs

cell = 3
a = cell_to_edges.ptrs[cell]
b = cell_to_edges.ptrs[cell+1]-1
edges = cell_to_edges.list[a:b]
#@show edges

# Generate a toy structured grid of the unit cube
# with 2x3x2 quadrilateral cells
grid = UGrid(domain=(0,1,0,1,0,1),partition=(2,3,2))

# Get the connectivity of the cells of the grid (i.e., for each cell, the ids of its vertices)
cell_to_vertices = connections(grid)

# The cell connectivity is represented by a `Connections` object
@assert isa(cell_to_vertices, Connections)

# `Connections` is a struct representing a vector of vectors
# in compressed form. It is fully described by a vector containing
# the flatted underlying data and another vector containing
# the start and end of each of the sub-vectors.
# That is:

cell_to_vertices_data = list(cell_to_vertices)
cell_to_vertices_ptrs = ptrs(cell_to_vertices)

cell = 2
ibeg = cell_to_vertices_ptrs[cell]
iend = cell_to_vertices_ptrs[cell+1]-1
vertices = cell_to_vertices_data[ibeg:iend] # Vertices of cell 2
@assert vertices == [2, 3, 5, 6, 14, 15, 17, 18]

# Get the coordinates of the mesh vertices
vertex_to_coords = coordinates(grid)
@assert isa(vertex_to_coords, Array{Float64,2})

# First axis for space dimensions, second one for number of vertices
@assert size(vertex_to_coords) == (3,36)

# Find cells around each vertex
vertex_to_cells = generate_dual_connections(cell_to_vertices)

# The returned data is a vector of vectors represented by a Connections object
@assert isa(vertex_to_cells, Connections)

# Generate a global numbering for the edges (1d objects) of the grid,
# and find which edges are on the boundary of each cell
n=1
cell_to_edges = generate_cell_to_faces(n,grid)

# The returned data again a Connections object
@assert isa(cell_to_edges, Connections)

# Idem for faces (2d objects)
n=2
cell_to_faces = generate_cell_to_faces(n,grid)
@assert isa(cell_to_faces, Connections)


text:
using JuMP
using Combinatorics

# Inspired from Joey Huchette's test in ConvexHull.jl
function nonfulldimensionaltest(lib::Polyhedra.Library)
    m = Model()
    @variable(m, x[1:3] ≥ 1)
    @constraints(m, begin
        x[1] == 2
        x[2] ≤ 2
    end)

    poly = polyhedron(m, lib)

    V = [2 1 1;
         2 2 1]
    R = [0 0 1]
    generator_fulltest(poly, V, R)
end


text:
struct Signature{P,N,D} end

Base.broadcastable(x::Signature) = Ref(x)

Signature(positive, negative=0, degenerate=0) = Signature{positive, negative, degenerate}()
Signature(string::AbstractString) = Signature(count.(["+", "-", "𝟎"], Ref(string))...)

positive(::Signature{P}) where {P} = P

negative(::Signature{P,N}) where {P,N} = N

degenerate(::Signature{P,N,D}) where {P,N,D} = D

dimension(::Signature{P,N,D}) where {P,N,D} = P + N + D

is_degenerate(sig::Signature) = degenerate(sig) ≠ 0

triplet(sig::Signature) = (positive(sig), negative(sig), degenerate(sig))

metric(::Signature{P,N,D}, ::Val{I}) where {P,N,D,I} = I <= P ? 1 : I <= P + N ? -1 : 0
metric(sig::Signature{P,N,D}, i::Val{I}, j::Val{I}) where {P,N,D,I} = metric(sig, i)
metric(::Signature, ::Val{I}, ::Val{J}) where {I,J} = 0

show(io::IO, sig::Signature) = print(io, sig == Ø ? "Ø" : "<" * join(["+", "-", "𝟎"] .^ [positive(sig), negative(sig), degenerate(sig)]) * ">")


text:
const Optional{T} = Union{Nothing, T}

struct RectSides{T}
    left::T
    right::T
    bottom::T
    top::T
end

abstract type Side end

struct Left <: Side end
struct Right <: Side end
struct Top <: Side end
struct Bottom <: Side end
# for protrusion content:
struct TopLeft <: Side end
struct TopRight <: Side end
struct BottomLeft <: Side end
struct BottomRight <: Side end

struct Inner <: Side end
struct Outer <: Side end

abstract type GridDir end
struct Col <: GridDir end
struct Row <: GridDir end

struct RowCols{T <: Union{Number, Vector{Float64}}}
    lefts::T
    rights::T
    tops::T
    bottoms::T
end

"""
    struct Span

Used to specify space that is occupied in a grid. Like 1:1|1:1 for the first square,
or 2:3|1:4 for a rect over the 2nd and 3rd row and the first four columns.
"""
struct Span
    rows::UnitRange{Int64}
    cols::UnitRange{Int64}
end

"""
    mutable struct GridContent{G, T}

Wraps content elements of a `GridLayout`. It keeps track of the `parent`, the `content` and its position in the grid via `span` and `side`.
"""
mutable struct GridContent{G, T} # G should be GridLayout but can't be used before definition
    parent::Optional{G}
    content::T
    span::Span
    side::Side
    needs_update::Observable{Bool}
    protrusions_handle::Optional{Function}
    reportedsize_handle::Optional{Function}
end

abstract type AlignMode end

"AlignMode that excludes the protrusions from the bounding box."
struct Inside <: AlignMode end

"AlignMode that includes the protrusions within the bounding box, plus paddings."
struct Outside <: AlignMode
    padding::RectSides{Float32}
end
Outside() = Outside(0f0)
Outside(padding::Real) = Outside(RectSides{Float32}(padding, padding, padding, padding))
Outside(left::Real, right::Real, bottom::Real, top::Real) =
    Outside(RectSides{Float32}(left, right, bottom, top))

"AlignMode that is Inside where padding is Nothing and Outside where it is Real."
struct Mixed <: AlignMode
    padding::RectSides{Union{Nothing, Float32}}
end
function Mixed(; left = nothing, right = nothing, bottom = nothing, top = nothing)
    paddings = map((left, right, bottom, top)) do side
        isnothing(side) ? side : Float32(side)
    end
    Mixed(RectSides{Union{Nothing, Float32}}(paddings...))
end

abstract type ContentSize end
abstract type GapSize <: ContentSize end

"""
    struct Auto <: ContentSize

If used as a `GridLayout`'s row / column size and `trydetermine == true`, signals to the `GridLayout` that the row / column should shrink to match the largest determinable element inside.
If no size of a content element can be determined, the remaining space is split between all `Auto` rows / columns according to their `ratio`.

If used as width / height of a layoutable element and `trydetermine == true`, the element's computed width / height will report the auto width / height if it can be determined.
This enables a parent `GridLayout` to adjust its column / rowsize to the element's width / height.
If `trydetermine == false`, the element's computed width / height will report `nothing` even if an auto width / height can be determined, which will prohibit a parent `GridLayout` from adjusting a row / column to the element's width / height.
This is useful to, e.g., prohibit a `GridLayout` from shrinking a column's width to the width of a super title, even though the title's width can be auto-determined.

The `ratio` is ignored if `Auto` is used as an element size.
"""
struct Auto <: ContentSize
    trydetermine::Bool # false for determinable size content that should be ignored
    ratio::Float64 # float ratio in case it's not determinable

    Auto(trydetermine::Bool = true, ratio::Real = 1.0) = new(trydetermine, ratio)
end
Auto(ratio::Real) = Auto(true, ratio)

struct Fixed <: GapSize
    x::Float64
end
struct Relative <: GapSize
    x::Float64
end
struct Aspect <: ContentSize
    index::Int
    ratio::Float64
end

"""
    mutable struct LayoutObservables{T, G}

`T` is the same type parameter of contained `GridContent`, `G` is `GridLayout` which is defined only after `LayoutObservables`.

A collection of `Observable`s and an optional `GridContent` that are needed to interface with the MakieLayout layouting system.

- `suggestedbbox::Observable{FRect2D}`: The bounding box that an element should place itself in. Depending on the element's `width` and `height` attributes, this is not necessarily equal to the computedbbox.
- `protrusions::Observable{RectSides{Float32}}`: The sizes of content "sticking out" of the main element into the `GridLayout` gaps.
- `reportedsize::Observable{NTuple{2, Optional{Float32}}}`: The width and height that the element computes for itself if possible (else `nothing`).
- `autosize::Observable{NTuple{2, Optional{Float32}}}`: The width and height that the element reports to its parent `GridLayout`. If the element doesn't want to cause the parent to adjust to its size, autosize can hide the reportedsize from it by being set to `nothing`.
- `computedbbox::Observable{FRect2D}`: The bounding box that the element computes for itself after it has received a suggestedbbox.
- `gridcontent::Optional{GridContent{G, T}}`: A reference of a `GridContent` if the element is currently placed in a `GridLayout`. This can be used to retrieve the parent layout, remove the element from it or change its position, and assign it to a different layout.
"""
mutable struct LayoutObservables{T, G} # G again GridLayout
    suggestedbbox::Observable{FRect2D}
    protrusions::Observable{RectSides{Float32}}
    reportedsize::Observable{NTuple{2, Optional{Float32}}}
    autosize::Observable{NTuple{2, Optional{Float32}}}
    computedbbox::Observable{FRect2D}
    gridcontent::Optional{GridContent{G, T}} # the connecting link to the gridlayout
end

mutable struct GridLayout
    content::Vector{GridContent}
    nrows::Int
    ncols::Int
    rowsizes::Vector{ContentSize}
    colsizes::Vector{ContentSize}
    addedrowgaps::Vector{GapSize}
    addedcolgaps::Vector{GapSize}
    alignmode::AlignMode
    equalprotrusiongaps::Tuple{Bool, Bool}
    needs_update::Observable{Bool}
    block_updates::Bool
    layoutobservables::LayoutObservables
    width::Observable
    height::Observable
    tellwidth::Observable
    tellheight::Observable
    halign::Observable
    valign::Observable
    default_rowgap::GapSize
    default_colgap::GapSize
    _update_func_handle::Optional{Function} # stores a reference to the result of on(obs)

    function GridLayout(
        content, nrows, ncols, rowsizes, colsizes,
        addedrowgaps, addedcolgaps, alignmode, equalprotrusiongaps, needs_update,
        layoutobservables, width, height, tellwidth, tellheight, halign, valign, default_rowgap, default_colgap)

        gl = new(content, nrows, ncols, rowsizes, colsizes,
            addedrowgaps, addedcolgaps, alignmode, equalprotrusiongaps,
            needs_update, false, layoutobservables, width, height, tellwidth, tellheight,
            halign, valign, default_rowgap, default_colgap, nothing)

        validategridlayout(gl)

        gl
    end
end

const Indexables = Union{UnitRange, Int, Colon}

struct GridPosition
    layout::GridLayout
    span::Span
    side::Side
end

struct GridLayoutSpec
    content::Vector{Pair{Tuple{Indexables, Indexables, Side}, Any}}
    kwargs::Dict{Symbol, Any}
end


text:
using Allnix: Point, norm, distance

@testset "throw" begin
@test_throws DomainError throw(DomainError())
@test_throws ErrorException error("I am an ErrorException")
# try
#     error("Error")
# catch e
#     println(e)
# end
end

@testset "isa" begin
@test isa(13.0, Float64)
@test isa(13, Int64)
end

@testset "types" begin
@test 3 == (1+2)::Int64
@test_throws TypeError 3 == (1.0 + 2.0)::Int64
end

@testset "Point" begin
a::Point{Float64} = Point{Float64}(1,1)
b::Point{Float64} = Point{Float64}(2,2)

@test sqrt(2) == norm(a)
@test sqrt(8) == norm(b)

@test sqrt(2) == distance(a,b)
end

@testset "Vector" begin
a = Vector{Float64}()
length = 101

@test Base.length(a) == 0
@test isempty(a)

resize!(a, length)

@test Base.length(a) == length

# - zero an array - #
a[:] = zero(Float64)

for x in a
    @test x == 0
end

a[1] = 1325
a[length] = 0.123
@test a[1] == 1325
@test a[length] == 0.123

# - another way to zero an array - #
fill!(a, 0.::Float64)

for x in a
    @test x == 0
end

# - yet another way - #

a[1] = 1325
a[length] = 0.123
@test a[1] == 1325
@test a[length] == 0.123

# TODO: broadcast!
a .= 0
for x in a
    @test x == 0
end

end


text:
"""
    ib_briowu(set, ps, vs, gas)
    ib_briowu(gam, mi, me, uspace)
    ib_briowu(gam, mi, me, uspace, vspace)

Initialize Brio-Wu MHD shock tube

"""
function ib_briowu(
    set::AbstractSetup,
    ps::AbstractPhysicalSpace,
    vs::Union{AbstractVelocitySpace,Nothing},
    gas::AbstractProperty,
)

    # upstream
    primL = zeros(5, 2)
    primL[1, 1] = 1.0 * gas.mi
    primL[2, 1] = 0.0
    primL[3, 1] = 0.0
    primL[4, 1] = 0.0
    primL[5, 1] = gas.mi / 1.0
    primL[1, 2] = 1.0 * gas.me
    primL[2, 2] = 0.0
    primL[3, 2] = 0.0
    primL[4, 2] = 0.0
    primL[5, 2] = gas.me / 1.0

    wL = mixture_prim_conserve(primL, gas.γ)

    EL = zeros(3)
    BL = zeros(3)
    BL[1] = 0.75
    BL[2] = 1.0
    lorenzL = zeros(3, 2)

    # downstream
    primR = zeros(5, 2)
    primR[1, 1] = 0.125 * gas.mi
    primR[2, 1] = 0.0
    primR[3, 1] = 0.0
    primR[4, 1] = 0.0
    primR[5, 1] = gas.mi * 1.25
    primR[1, 2] = 0.125 * gas.me
    primR[2, 2] = 0.0
    primR[3, 2] = 0.0
    primR[4, 2] = 0.0
    primR[5, 2] = gas.me * 1.25

    wR = mixture_prim_conserve(primR, gas.γ)
    ER = zeros(3)
    BR = zeros(3)
    BR[1] = 0.75
    BR[2] = -1.0
    lorenzR = zeros(3, 2)

    fw = function (x)
        if x <= (ps.x0 + ps.x1) / 2
            return wL
        else
            return wR
        end
    end
    fE = function (x)
        if x <= (ps.x0 + ps.x1) / 2
            return EL
        else
            return ER
        end
    end
    fB = function (x)
        if x <= (ps.x0 + ps.x1) / 2
            return BL
        else
            return BR
        end
    end
    fL = function (x)
        if x <= (ps.x0 + ps.x1) / 2
            return lorenzL
        else
            return lorenzR
        end
    end

    bc = function (x)
        if x <= (ps.x0 + ps.x1) / 2
            return primL
        else
            return primR
        end
    end

    if set.space[3:end] == "4f1v"
        h0L = mixture_maxwellian(vs.u, primL)
        h1L = similar(h0L)
        h2L = similar(h0L)
        h3L = similar(h0L)
        for j in axes(h0L, 2)
            h1L[:, j] .= primL[3, j] .* h0L[:, j]
            h2L[:, j] .= primL[4, j] .* h0L[:, j]
            h3L[:, j] .=
                (primL[3, j]^2 + primL[4, j]^2 + 2.0 / (2.0 * primL[end, j])) .* h0L[:, j]
        end

        h0R = mixture_maxwellian(vs.u, primR)
        h1R = similar(h0R)
        h2R = similar(h0R)
        h3R = similar(h0R)
        for j in axes(h0L, 2)
            h1R[:, j] .= primR[3, j] .* h0R[:, j]
            h2R[:, j] .= primR[4, j] .* h0R[:, j]
            h3R[:, j] .=
                (primR[3, j]^2 + primR[4, j]^2 + 2.0 / (2.0 * primR[end, j])) .* h0R[:, j]
        end

        ff = function (x)
            if x <= (ps.x0 + ps.x1) / 2
                return h0L, h1L, h2L, h3L
            else
                return h0R, h1R, h2R, h3R
            end
        end

        return fw, ff, fE, fB, fL, bc
    elseif set.space[3:end] == "3f2v"
        h0L = mixture_maxwellian(vs.u, vs.v, primL)
        h1L = similar(h0L)
        h2L = similar(h0L)
        for j in axes(h0L, 3)
            h1L[:, :, j] .= primL[4, j] .* h0L[:, :, j]
            h2L[:, :, j] .= (primL[4, j]^2 + 1.0 / (2.0 * primL[end, j])) .* h0L[:, :, j]
        end

        h0R = mixture_maxwellian(vs.u, vs.v, primR)
        h1R = similar(h0R)
        h2R = similar(h0R)
        for j in axes(h0R, 3)
            h1R[:, :, j] .= primR[4, j] .* h0R[:, :, j]
            h2R[:, :, j] .= (primR[4, j]^2 + 1.0 / (2.0 * primR[end, j])) .* h0R[:, :, j]
        end

        ff = function (x)
            if x <= (ps.x0 + ps.x1) / 2
                return h0L, h1L, h2L
            else
                return h0R, h1R, h2R
            end
        end

        return fw, ff, fE, fB, fL, bc
    end

    return nothing

end

function ib_briowu(gam, mi, me, uspace::AbstractMatrix{T}) where {T<:AbstractFloat}

    # upstream
    primL = zeros(5, 2)
    primL[1, 1] = 1.0 * mi
    primL[2, 1] = 0.0
    primL[3, 1] = 0.0
    primL[4, 1] = 0.0
    primL[5, 1] = mi / 1.0
    primL[1, 2] = 1.0 * me
    primL[2, 2] = 0.0
    primL[3, 2] = 0.0
    primL[4, 2] = 0.0
    primL[5, 2] = me / 1.0

    wL = mixture_prim_conserve(primL, gam)
    h0L = mixture_maxwellian(uspace, primL)

    h1L = similar(h0L)
    h2L = similar(h0L)
    h3L = similar(h0L)
    for j in axes(h0L, 2)
        h1L[:, j] .= primL[3, j] .* h0L[:, j]
        h2L[:, j] .= primL[4, j] .* h0L[:, j]
        h3L[:, j] .=
            (primL[3, j]^2 + primL[4, j]^2 + 2.0 / (2.0 * primL[end, j])) .* h0L[:, j]
    end

    EL = zeros(3)
    BL = zeros(3)
    BL[1] = 0.75
    BL[2] = 1.0

    # downstream
    primR = zeros(5, 2)
    primR[1, 1] = 0.125 * mi
    primR[2, 1] = 0.0
    primR[3, 1] = 0.0
    primR[4, 1] = 0.0
    primR[5, 1] = mi * 1.25
    primR[1, 2] = 0.125 * me
    primR[2, 2] = 0.0
    primR[3, 2] = 0.0
    primR[4, 2] = 0.0
    primR[5, 2] = me * 1.25

    wR = mixture_prim_conserve(primR, gam)
    h0R = mixture_maxwellian(uspace, primR)

    h1R = similar(h0R)
    h2R = similar(h0R)
    h3R = similar(h0R)
    for j in axes(h0L, 2)
        h1R[:, j] .= primR[3, j] .* h0R[:, j]
        h2R[:, j] .= primR[4, j] .* h0R[:, j]
        h3R[:, j] .=
            (primR[3, j]^2 + primR[4, j]^2 + 2.0 / (2.0 * primR[end, j])) .* h0R[:, j]
    end

    ER = zeros(3)
    BR = zeros(3)
    BR[1] = 0.75
    BR[2] = -1.0

    lorenzL = zeros(3, 2)
    lorenzR = zeros(3, 2)

    return wL,
    primL,
    h0L,
    h1L,
    h2L,
    h3L,
    EL,
    BL,
    lorenzL,
    wR,
    primR,
    h0R,
    h1R,
    h2R,
    h3R,
    ER,
    BR,
    lorenzR

end

function ib_briowu(
    gam,
    mi,
    me,
    uspace::T,
    vspace::T,
) where {T<:AbstractArray{<:AbstractFloat,3}}

    # upstream
    primL = zeros(5, 2)
    primL[1, 1] = 1.0 * mi
    primL[2, 1] = 0.0
    primL[3, 1] = 0.0
    primL[4, 1] = 0.0
    primL[5, 1] = mi / 1.0
    primL[1, 2] = 1.0 * me
    primL[2, 2] = 0.0
    primL[3, 2] = 0.0
    primL[4, 2] = 0.0
    primL[5, 2] = me / 1.0

    wL = mixture_prim_conserve(primL, gam)
    h0L = mixture_maxwellian(uspace, vspace, primL)

    h1L = similar(h0L)
    h2L = similar(h0L)
    for j in axes(h0L, 3)
        h1L[:, :, j] .= primL[4, j] .* h0L[:, :, j]
        h2L[:, :, j] .= (primL[4, j]^2 + 1.0 / (2.0 * primL[end, j])) .* h0L[:, :, j]
    end

    EL = zeros(3)
    BL = zeros(3)
    BL[1] = 0.75
    BL[2] = 1.0

    # downstream
    primR = zeros(5, 2)
    primR[1, 1] = 0.125 * mi
    primR[2, 1] = 0.0
    primR[3, 1] = 0.0
    primR[4, 1] = 0.0
    primR[5, 1] = mi * 1.25
    primR[1, 2] = 0.125 * me
    primR[2, 2] = 0.0
    primR[3, 2] = 0.0
    primR[4, 2] = 0.0
    primR[5, 2] = me * 1.25

    wR = mixture_prim_conserve(primR, gam)
    h0R = mixture_maxwellian(uspace, vspace, primR)

    h1R = similar(h0R)
    h2R = similar(h0R)
    for j in axes(h0R, 3)
        h1R[:, :, j] .= primR[4, j] .* h0R[:, :, j]
        h2R[:, :, j] .= (primR[4, j]^2 + 1.0 / (2.0 * primR[end, j])) .* h0R[:, :, j]
    end

    ER = zeros(3)
    BR = zeros(3)
    BR[1] = 0.75
    BR[2] = -1.0

    lorenzL = zeros(3, 2)
    lorenzR = zeros(3, 2)

    return wL,
    primL,
    h0L,
    h1L,
    h2L,
    EL,
    BL,
    lorenzL,
    wR,
    primR,
    h0R,
    h1R,
    h2R,
    ER,
    BR,
    lorenzR

end


text:
using Test
using Pretend
Pretend.activate()

# Mocking.jl example
@mockable open(f,io) = Base.open(f,io)
@mockable read(fp, n) = Base.read(fp, n)

function randdev(n::Integer)
    open("/dev/urandom") do fp
        reverse(read(fp, n))
    end
end

# patch `read` to test properly
# patch `open` to avoid opening any file
@test apply(read => (fp,n) -> [UInt8(i) for i in 1:n],
            open => (f,io) -> f(io)) do
    randdev(5)
end == [0x05, 0x04, 0x03, 0x02, 0x01]

# outside test it's a no-op
@test randdev(5) != [0x05, 0x04, 0x03, 0x02, 0x01]


text:
using QuantizedStateSystems

@model function my_model(t, x)
  dx = t*x
end


text:
import Base.*, Base.+, Base.-, Base./, Base.show, Base.!=, Base.==, Base.<=, Base.<, Base.>, Base.>=, Base.divrem

const z0 = "0"
const z1 = "1"
const flipordered = (z1 < z0)

mutable struct Z s::String end
Z() = Z(z0)
Z(z::Z) = Z(z.s)

pairlen(x::Z, y::Z) = max(length(x.s), length(y.s))
tolen(x::Z, n::Int) = (s = x.s; while length(s) < n s = z0 * s end; s)

<(x::Z, y::Z) = (l = pairlen(x, y); flipordered ? tolen(x, l) > tolen(y, l) : tolen(x, l) < tolen(y, l))
>(x::Z, y::Z) = (l = pairlen(x, y); flipordered ? tolen(x, l) < tolen(y, l) : tolen(x, l) > tolen(y, l))
==(x::Z, y::Z) = (l = pairlen(x, y); tolen(x, l) == tolen(y, l))
<=(x::Z, y::Z) = (l = pairlen(x, y); flipordered ? tolen(x, l) >= tolen(y, l) : tolen(x, l) <= tolen(y, l))
>=(x::Z, y::Z) = (l = pairlen(x, y); flipordered ? tolen(x, l) <= tolen(y, l) : tolen(x, l) >= tolen(y, l))
!=(x::Z, y::Z) = (l = pairlen(x, y); tolen(x, l) != tolen(y, l))

function tocanonical(z::Z)
    while occursin(z0 * z1 * z1, z.s)
        z.s = replace(z.s, z0 * z1 * z1 => z1 * z0 * z0)
    end
    len = length(z.s)
    if len > 1 && z.s[1:2] == z1 * z1
        z.s = z1 * z0 * z0 * ((len > 2) ? z.s[3:end] : "")
    end
    while (len = length(z.s)) > 1 && string(z.s[1]) == z0
        if len == 2
            if z.s == z0 * z0
                z.s = z0
            elseif z.s == z0 * z1
                z.s = z1
            end
        else
            z.s = z.s[2:end]
        end
    end
    z
end

function inc(z)
    if z.s[end] == z0[1]
        z.s = z.s[1:end-1] * z1[1]
    elseif z.s[end] == z1[1]
        if length(z.s) > 1
            if z.s[end-1:end] == z0 * z1
                z.s = z.s[1:end-2] * z1 * z0
            end
        else
            z.s = z1 * z0
        end
    end
    tocanonical(z)
end

function dec(z)
    if z.s[end] == z1[1]
        z.s = z.s[1:end-1] * z0
    else
        if (m = match(Regex(z1 * z0 * '+' * '$'), z.s)) != nothing
            len = length(m.match)
            if iseven(len)
                z.s = z.s[1:end-len] * (z0 * z1) ^ div(len, 2)
            else
                z.s = z.s[1:end-len] * (z0 * z1) ^ div(len, 2) * z0
            end
        end
    end
    tocanonical(z)
    z
end

function +(x::Z, y::Z)
    a = Z(x.s)
    b = Z(y.s)
    while b.s != z0
        inc(a)
        dec(b)
    end
    a
end

function -(x::Z, y::Z)
    a = Z(x.s)
    b = Z(y.s)
    while b.s != z0
        dec(a)
        dec(b)
    end
    a
end

function *(x::Z, y::Z)
    if (x.s == z0) || (y.s == z0)
        return Z(z0)
    elseif x.s == z1
        return Z(y.s)
    elseif y.s == z1
        return Z(x.s)
    end
    a = Z(x.s)
    b = Z(z1)
    while b != y
        c = Z(z0)
        while c != x
            inc(a)
            inc(c)
        end
        inc(b)
    end
    a
end

function divrem(x::Z, y::Z)
    if y.s == z0
        throw("Zeckendorf division by 0")
    elseif (y.s == z1) || (x.s == z0)
        return Z(x.s)
    end
    a = Z(x.s)
    b = Z(y.s)
    c = Z(z0)
    while a > b
        a = a - b
        inc(c)
    end
    tocanonical(c), tocanonical(a)
end

function /(x::Z, y::Z)
    a, _ = divrem(x, y)
    a
end

show(io::IO, z::Z) = show(io, parse(BigInt, tocanonical(z).s))

function zeckendorftest()
    a = Z("10")
    b = Z("1001")
    c = Z("1000")
    d = Z("10101")

    println("Addition:")
    x = a
    println(x += a)
    println(x += a)
    println(x += b)
    println(x += c)
    println(x += d)

    println("\nSubtraction:")
    x = Z("1000")
    println(x - Z("101"))
    x = Z("10101010")
    println(x - Z("1010101"))

    println("\nMultiplication:")
    x = Z("1001")
    y = Z("101")
    println(x * y)
    println(Z("101010") * y)

    println("\nDivision:")
    x = Z("1000101")
    y = Z("101")
    println(x / y)
    println(divrem(x, y))
end

zeckendorftest()

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode