text:
module SYM

using BenchmarkTools: push!
using BenchmarkTools

import Base:+,-,*,/,^,sin,cos,tan,exp # methods for ordinary numbers are needed, so we have to import explicitly (since we have changed their methods)
import Base:convert, show

export Sym, show
export @def, def_from_vec
export sub!, lambdify, lambdify_fast

#####################################################################
########## Types, Promotion Rules, and Convert Start Here ###########
#####################################################################
"""
An instance of Symbolic variable can be realized by
	`Sym{T}(expr)`, where T is the SubType of the Real domain.
Here I confine the variable in real just for convenience of the promotion rule with complex numbers since Complex{T}(x) is defined only for x of the
	subtype of Real, while Complex{Sym} is illegal.
"""
abstract type Symbolic{T} end # <: Number

struct Sym{T<:Number} <: Symbolic{T}
	expr # of Type{Symbol} or Type{Expr}
end

"Claim realization of the instance for other types"
Sym{T}(x::String) where {T<:Number} = Sym{T}(Meta.parse(x))
Sym{T}(x::Sym) where {T<:Number} = Sym{T}(x.expr) # change the symtype, or keep unchanged for repeated invokation
Sym(x) = Sym{Real}(x) # add a method for default type of Real
Sym(list::Vector) = [Sym{Real}(i) for i in list] # add a method for generating a list of symbols
#= for example,
	x = Sym{Real}(:x) # creat one instance of Real Symbolic variable
	y = Sym{Real}("y")
	z = Sym(x)
	x = Sym{Real}(gensym())	# creat one instance for an anonymous variable
=#

"Define a new type-check function by adding methods"
function symtype end
symtype(::Sym{T}) where {T} = T
symtype(x) = typeof(x) # for other known Types of x like Number, Symbol, and Expr etc.

"""
Define promote rules for Type{Sym} and Type{Number}.
This is necessary since basic binary operations always act in the way of, for example, +(x::Number, y::Number) = +(promote(x,y)...)
See https://docs.julialang.org/en/v1/manual/conversion-and-promotion/#Promotion
"""
Base.promote_rule(x::Type{Sym{T}}, y::Type{S}) where {T, S<:Real} = Sym{Union{T, S}}
# for example, promote_rule(Sym{Float64}, Int64) = Sym{Float64}

"""
Add methods for basic binary operations +,-,*,/,^
Here we distinguish real and complex numbers
"""
+(x::Sym{T}, y::Number) where {T} = Sym{T}(eval(:(:($($x.expr)+$$y))))
+(y::Number, x::Sym{T}) where {T} = Sym{T}(eval(:(:($($x.expr)+$$y))))
+(x::Sym{T}, y::Sym{S}) where {T,S} = Sym{Union{T,S}}(eval(:(:($($x.expr)+$($y.expr)))))
-(x::Sym{T}, y::Number) where {T} = Sym{T}(eval(:(:($($x.expr)-$$y))))
-(y::Number, x::Sym{T}) where {T} = Sym{T}(eval(:(:($($x.expr)-$$y))))
-(x::Sym{T}, y::Sym{S}) where {T,S} = Sym{Union{T,S}}(eval(:(:($($x.expr)-$($y.expr)))))
*(x::Sym{T}, y::Number) where {T} = Sym{T}(eval(:(:($($x.expr)*$$y))))
*(y::Number, x::Sym{T}) where {T} = Sym{T}(eval(:(:($($x.expr)*$$y))))
*(x::Sym{T}, y::Sym{S}) where {T,S} = Sym{Union{T,S}}(eval(:(:($($x.expr)*$($y.expr)))))
/(x::Sym{T}, y::Number) where {T} = Sym{T}(eval(:(:($($x.expr)/$$y))))
/(y::Number, x::Sym{T}) where {T} = Sym{T}(eval(:(:($$y/$($x.expr)))))
/(x::Sym{T}, y::Sym{S}) where {T,S} = Sym{Union{T,S}}(eval(:(:($($x.expr)/$($y.expr)))))
^(x::Sym{T}, y::Number) where {T} = Sym{T}(eval(:(:($($x.expr)^$$y))))
^(y::Number, x::Sym{T}) where {T} = Sym{T}(eval(:(:($$y^$($x.expr)))))
^(x::Sym{T}, y::Sym{S}) where {T,S} = Sym{Union{T,S}}(eval(:(:($($x.expr)^$($y.expr)))))

# define the zero unit for abelian groups
Base.one(s::Sym) = one(symtype(s))
Base.zero(s::Sym) = zero(symtype(s))

"""
Add methods for basic single-variable mathematical functions exp, sin, cos, tan (to be continue...)
"""
exp(x::Sym{T}) where {T} = Sym{T}(eval(:(:(exp($($x.expr))))))
exp(x::Array{Sym{T}}) where {T} = Sym{T}(:(exp($x)))
sin(x::Sym{T}) where {T} = Sym{T}(eval(:(:(sin($($x.expr))))))
sin(x::Array{Sym{T}}) where {T} = Sym{T}(:(sin($x)))
cos(x::Sym{T}) where {T} = Sym{T}(eval(:(:(cos($($x.expr))))))
cos(x::Array{Sym{T}}) where {T} = Sym{T}(:(cos($x)))
tan(x::Sym{T}) where {T} = Sym{T}(eval(:(:(tan($($x.expr))))))
tan(x::Array{Sym{T}}) where {T} = Sym{T}(:(tan($x)))

"""
Convert methods
	motivated by https://github.com/symengine/SymEngine.jl/blob/master/src/subs.jl
"""
function convert(::Type{Expr}, x::Sym{T}) where {T}
	return x.expr
end

function convert(::Type{Expr}, v::Vector{Sym{T}}) where {T}# slow for large array
    Expr(:vcat, map(x->x.expr, v)...) # map .expr to the vector v. Such Vectorization treatment is much faster.
end

function convert(::Type{Expr}, m::Matrix{Sym{T}}) where {T} # slow for large array
    col_args = []
    for i = 1:size(m, 1)
        row_args = []
        for j = 1:size(m, 2)
            push!(row_args, m[i,j].expr)
        end
        row = Expr(:hcat, row_args...)
        push!(col_args, row)
    end
    Expr(:vcat, col_args...)
end

"Change the show methods for Type{Sym} (default looks like Sym{T}(:x); here we now display its Expr representation"
function Base.show(io::IO, x::Sym)
	print(io, "$(x.expr)")
end

#####################################################################
######################### Macros Start Here #########################
#####################################################################
"""
Claim (real) symbolic variables (with the Type{Sym{Real}}) by its Symbol representations
"""
macro def(x...) # accept multiple arguments
	q = Expr(:block)
	for s in x
		if s isa String
			s = Meta.parse(s) # convert to Symbol or Expr
			push!(q.args, Expr(:(=), esc(s), Expr(:call, :Sym, Expr(:quote, s)))) # Expr for claiming symbolic variables
		elseif s isa Symbol
			push!(q.args, Expr(:(=), esc(s), Expr(:call, :Sym, Expr(:quote, s)))) # Expr for claiming symbolic variables
		else
			error("Symbolic variables cannot be constructed: Not implemented!")
		end
	end
	push!(q.args, Expr(:vect, map(esc, x)...)) # the last argument of q is a vector of all variables.
	return q # automatically evaluate the Exprs so that all variables are claimed globally. And the last argument be evaluated and displayed.
end
#=
"Claim (real) symbolic variables (with the Type{Sym{Real}}) by its Symbol representations"
function def_from_vec(x::Vector) # accept multiple arguments
	q = []
	for s in x
		if isa(s, Symbol)
			eval(Expr(:(=), s, Expr(:call, :Sym, Expr(:quote, s)))) # claim symbolic variables
			push!(q, s)
		elseif isa(s, String)
			s = Meta.parse(s)
			eval(Expr(:(=), s, Expr(:call, :Sym, Expr(:quote, s)))) # claim symbolic variables
			push!(q, s)
		end
	end
	return q
end=#

#####################################################################
########### Substitution and Lambdification Start Here ##############
#####################################################################
const binary_op_list = [:+, :-, :*, :/, :^]

function AST_iterate_search_replacement(ex, rules::Tuple) # Slow serial tree search, wait to be optimized
	if ex isa Expr
		if ex.head == :call
			for (i,s) in enumerate(ex.args)
				if s isa Symbol && !(s in binary_op_list)
					for r in rules
						if s == r.first.expr
							if r.second isa Sym
								ex.args[i] = r.second.expr
							elseif r.second isa Union{Number, Symbol}
								ex.args[i] = r.second
							else
								error("Illegal replacement for the second arguments!")
							end
						end
					end
				elseif s isa Expr
					AST_iterate_search_replacement(s, rules)
				end
			end
		else
			error("Illegal expression: not implemented!")
		end
	elseif ex isa Symbol
		for r in rules
			if ex == r.first.expr
				if r.second isa Sym
					ex = r.second.expr
				elseif r.second isa Union{Number, Symbol}
					ex = r.second
				else
					error("Illegal replacement for the second arguments!")
				end
			end
		end
	end
	return ex
end

function sub!(sym::Sym, rules::Tuple) # motivated by the discussion of https://discourse.julialang.org/t/substitute-symbols-in-expression/23099
	ex = sym.expr
	ex = AST_iterate_search_replacement(ex, rules)
	return ex # fortunately, `eval(x::Sym)=x`.
end

function sub!(ex, rules::Tuple) # ex isa Symbol or Expr is OK
	return AST_iterate_search_replacement(ex, rules)
end

function sub!(m::Array, rules::Tuple) # add the method for substituting Array{Sym} of any dimension
	return sub!.(m, Ref(rules))
end

"""
Lambdify a muti-variable function with arbitrary vector of Symbols
	Symbolic lambdification (with variables not given enough, or extra variables given) is also supported!
"""
function lambdify(x::Sym, vars::Vector{Symbol})
	return eval(Expr(:function, Expr(:call, gensym(), vars...), x.expr))
	# only Symbol can be legal function arguments, so we have to `convert(Expr, x)`, which is simply done by `x.expr`.
end

function lambdify(m::Array, vars::Vector{Symbol}) # add a method for lambdifying Vector{Sym} or Matrix{Sym}
	return eval(Expr(:function, Expr(:call, gensym(), vars...), convert(Expr, m))) # such convert is slow for large arrays. Wait to be optimized.
end

"""
Fast lambdification (with less memory cost) by direct substitution
	x accept both Type{Sym} and Type{Array{Sym}} of any dimensions.
"""
function lambdify_fast(x, vars::Vector{Sym{T}}) where {T}
	function func(args::Number...)
		eval.(sub!(x, Tuple(Pair.(vars, args))))
		# Vectorized evaluation
	end
	return func
end

#####################################################################
########################## Tests Start Here #########################
#####################################################################

end


text:
using FractalTools

A1 = reshape([1/2],1,1)
b1 = [0]
A2 = reshape([1/2],1,1)
b2 = [1/2]

w1 = Transformation(A1,b1)
w2 = Transformation(A2,b2)

w = [w1, w2]
ifs = IFS(w)

generator = randalg_sequential_generator(ifs.ws, ifs.probs)

take!(generator)

ext: jl


text:
export DenseSplit

struct DenseSplit{VT,MT,C}
    Wr::MT
    Wc::MT
    b::VT
    σ::C
end
functor(d::DenseSplit) = (Wr=d.Wr, Wc=d.Wc, b=d.b), t -> DenseSplit(t.Wr,t.Wc,t.b,d.σ)
(l::DenseSplit)((σr, σc)) = l.σ.(l.Wr*σr .+ l.Wc*σc .+ l.b)

DenseSplit(in::Integer, args...;kwargs...) =
    DenseSplit(Complex{STD_REAL_PREC}, in, args...;kwargs...)
function DenseSplit(T::Type, in::Integer, out::Integer, σ = identity;
               initW = glorot_uniform, initb = glorot_uniform)
  return DenseSplit(initW(T, out, in), initW(T, out, in), initb(T, out), σ)
end

function Base.show(io::IO, l::DenseSplit)
  print(io, "DenseSplit(", size(l.Wr, 2), ", ", size(l.Wr, 1))
  l.σ == identity || print(io, ", ", l.σ)
  print(io, ")")
end

struct DenseSplitCache{Ta,Tb,Tc,Td,Te}
    σr::Tc
    σc::Tc
    out::Tb

    δℒℒ::Td
    δℒr::Te
    δℒc::Te

    θ::Ta
    θ_tmp::Ta
    out2::Tb
    valid::Bool
end

cache(l::DenseSplit{Ta,Tb}, arr_T, in_T, in_sz) where {Ta,Tb} =
    DenseSplitCache(similar(l.Wr, size(l.Wr,2)),
               similar(l.Wr, size(l.Wr,2)),
               similar(l.b),

               similar(l.Wr, size(l.Wr,1)),
               similar(l.Wr, 1, size(l.Wr, 2)),
               similar(l.Wr, 1, size(l.Wr, 2)),

               similar(l.b),
               similar(l.b),
               similar(l.b),
               false)

function layer_out_type_size(l::DenseSplit, in_T ,in_sz)
    T1     = promote_type(in_T, eltype(l.Wr))
    out_T  = promote_type(T1, eltype(l.b))
    out_sz = size(l.b)
    return out_T, out_sz
end

function (l::DenseSplit)(c::DenseSplitCache, (xr, xc))
    # The preallocated caches
    out  = c.out
    θ = c.θ
    θ₂ = c.θ_tmp

    # Store the input to this layer for the backpropagation
    σr = copyto!(c.σr, xr)
    σc = copyto!(c.σc, xc)

    #θ .= net.b .+ net.W * x
    mul!(θ, l.Wr, σr)
    mul!(θ₂, l.Wc, σc)
    θ .+= l.b .+ θ₂

    # Apply the nonlinear function
    out  .= l.σ.(θ)
    return out
end

function backprop(∇, l::DenseSplit, c::DenseSplitCache, δℒ)
    # The preallocated caches
    θ = c.θ
    δℒℒ = c.δℒℒ

    # Compute the actual sensitivity
    copyto!(δℒℒ, δℒ)
    δℒℒ .*= fwd_der.(l.σ, θ)

    ∇.Wr .= δℒℒ.*transpose(c.σr)
    ∇.Wc .= δℒℒ.*transpose(c.σc)
    ∇.b .= δℒℒ

    # TODO alloc of transpose
    δℒr = mul!(c.δℒr, transpose(δℒℒ), l.Wr)
    δℒc = mul!(c.δℒc, transpose(δℒℒ), l.Wc)

    return (δℒr, δℒc)
end


text:
using BinaryBuilder, Pkg

name = "SSGraphBLAS"
version = v"7.0.3"

# Collection of sources required to build SuiteSparse:GraphBLAS
sources = [
    GitSource("https://github.com/DrTimothyAldenDavis/GraphBLAS.git",
        "4bb44bb93a30b54c7039d1c6f5d87e3c85ec45a5")
]

# Bash recipe for building across all platforms
script = raw"""
# Compile GraphBLAS
cd $WORKSPACE/srcdir/GraphBLAS

if [[ "$target" == *-mingw* ]]; then
    CMAKE_OPTIONS="-DGBNCPUFEAT=1"
fi
make -j${nproc} CMAKE_OPTIONS="${CMAKE_OPTIONS} -DCMAKE_INSTALL_PREFIX=${prefix} -DCMAKE_BUILD_TYPE=Release -DCMAKE_TOOLCHAIN_FILE=${CMAKE_TARGET_TOOLCHAIN}"
make install
install_license LICENSE
if [[ ! -f "${libdir}/libgraphblas.${dlext}" ]]; then
    # For mysterious reasons, the shared library is not installed
    # when building for Windows
    mkdir -p "${libdir}"
    cp "build/libgraphblas.${dlext}" "${libdir}"
fi
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = supported_platforms(;experimental=true)

# The products that we will ensure are always built
products = [
    LibraryProduct("libgraphblas", :libgraphblas),
]

# Dependencies that must be installed before this package can be built
dependencies = [
    Dependency("CompilerSupportLibraries_jll"),
]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies;
               preferred_gcc_version=v"9", julia_compat="1.6")

max_stars_repo_name: MichelJuillard/Yggdrasil

max_issues_repo_name: MichelJuillard/Yggdrasil

max_forks_repo_name: MichelJuillard/Yggdrasil


text:
import DelimitedFiles
import Plots

function populate(fname::String, n_days::Int; n_init=-1)
  internal_timers = [n_init]
  if n_init == -1
    internal_timers = vec(DelimitedFiles.readdlm(fname, ',', Int))
  end

  n_size = Vector{Int}()
  for n in 1:1:n_days
    internal_timers .-= 1

    # Nomenclature misleading here. At 'zero' time, the fish continues to live
    # on and only *the following day* it breeds a new fish.
    zero_timers = internal_timers .== -1
    n_zeros = sum(zero_timers)
    internal_timers[zero_timers] .= 6
    internal_timers = vcat(internal_timers, zeros(Int, n_zeros) .+ 8)
    push!(n_size, size(internal_timers, 1))
  end
  n_fish = size(internal_timers, 1)
  println("$n_fish fishs after $n_days days.")
end

function more_efficient(fname, n_days)
  initial_timers = vec(DelimitedFiles.readdlm(fname, ',', Int))
  n_timers = vec(zeros(Int, 9))

  for initial_timer in initial_timers
    n_timers[initial_timer+1] += 1
  end

  for day in 1:1:n_days
    n_new_fishs = n_timers[1]
    n_timers[1:1:end-1] = n_timers[2:1:end]
    n_timers[end] = n_new_fishs
    n_timers[7] += n_new_fishs
  end
  n_fish = sum(n_timers)
  println("$n_fish fishs after $n_days days.")
end

populate("test_input.txt", 80)
populate("input.txt", 80)

more_efficient("test_input.txt", 18)
more_efficient("test_input.txt", 80)
more_efficient("test_input.txt", 256)
more_efficient("input.txt", 256)


text:
# MIT License
#
# Copyright (c) 2018 Martin Biel
#
# Permission is hereby granted, free of charge, to any person obtaining a copy
# of this software and associated documentation files (the "Software"), to deal
# in the Software without restriction, including without limitation the rights
# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
# copies of the Software, and to permit persons to whom the Software is
# furnished to do so, subject to the following conditions:
#
# The above copyright notice and this permission notice shall be included in all
# copies or substantial portions of the Software.
#
# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
# SOFTWARE.

# Common
# ------------------------------------------------------------
function MOI.get(regularizer::AbstractRegularizer, param::RawRegularizationParameter)
    name = Symbol(param.name)
    if !(name in fieldnames(typeof(regularizer.parameters)))
        error("Unrecognized parameter name: $(name) for regularizer $(typeof(regularizer)).")
    end
    return getfield(regularizer.parameters, name)
end

function MOI.set(regularizer::AbstractRegularizer, param::RawRegularizationParameter, value)
    name = Symbol(param.name)
    if !(name in fieldnames(typeof(regularizer.parameters)))
        error("Unrecognized parameter name: $(name) for regularizer $(typeof(regularizer)).")
    end
    setfield!(regularizer.parameters, name, value)
    return nothing
end

function add_projection_targets!(regularization::AbstractRegularization, model::MOI.AbstractOptimizer)
    ξ = regularization.ξ
    for i in eachindex(ξ)
        name = add_subscript(:ξ, i)
        set = SingleDecisionSet(1, ξ[i], NoSpecifiedConstraint(), false)
        var_index, _ = MOI.add_constrained_variable(model, set)
        set_decision!(regularization.decisions, var_index, ξ[i])
        MOI.set(model, MOI.VariableName(), var_index, name)
        regularization.projection_targets[i] = var_index
    end
    return nothing
end

function decision(::AbstractLShaped, regularization::AbstractRegularization)
    return map(regularization.ξ) do ξᵢ
        ξᵢ.value
    end
end

function objective_value(::AbstractLShaped, regularization::AbstractRegularization)
    return regularization.data.Q̃
end

function gap(lshaped::AbstractLShaped, regularization::AbstractRegularization)
    @unpack θ = lshaped.data
    @unpack Q̃ = regularization.data
    return abs(θ-Q̃)/(abs(Q̃)+1e-10)
end

function process_cut!(lshaped::AbstractLShaped, cut::AbstractHyperPlane, ::AbstractRegularization)
    return nothing
end


text:
using Flux, Statistics
using Flux: onecold, crossentropy

include("myplots.jl")
include("model-commons.jl")

epochs = 300
eta = 0.001

# Train the ANN model
function train_ann(d::DataSet; hidden=750, epochs=epochs, eta=eta, plotloss=true)
  @info("Constructing ANN model...")
  m = Chain(
    Dense(size(d.X_train[:,1],1), hidden, relu),
    Dense(hidden, size(d.Y_train[:,1],1)),
    softmax)

  loss(x, y) = crossentropy(m(x), y)
  accuracy(x, y) = mean(onecold(m(x)) .== onecold(y))

  dataset = [(d.X_train, d.Y_train)]
  evalcb = () -> @show(loss(d.X_train, d.Y_train))
  opt = ADAM(eta)

  # Check for NaN values so we know if we have bad data before training
  # @info("Checking for NaN values")
  # for i in 1:size(X_train,2)
  #   if (isnan(loss(X_train[:,i], Y_train[:,i])))
  #     @info("$i is NaN")
  #   end
  # end
  # @info("Done.")

  @info("Training model...")

  # Train for a specified number of epochs
  if plotloss
    plt = create_plot(0, loss(d.X_train, d.Y_train).data)
  end
  for step in 1:epochs
    @info("Epoch $step")
    Flux.train!(loss, params(m), dataset, opt)
    error = loss(d.X_train, d.Y_train).data
    @show(error)
    @show(accuracy(d.X_test, d.Y_test))
    if plotloss
      update_plot(plt, step, error)
    end
  end

  @info("Done")

  return m
end

# 250 epochs, image size 150x100, eta 0.001, model = 15000 relu | 750 sigmoid | 10 softmax => 77.6% accuracy


text:
using PowerSystems
using JuMP

base_dir = string(dirname(dirname(pathof(PowerSystems))))
println(joinpath(base_dir,"data/data_5bus.jl"))
include(joinpath(base_dir,"data/data_5bus.jl"))

battery = GenericBattery(name = "Bat",
                status = true,
                energy = 10.0,
                activepower = 10.0,
                capacity = (min = 0.0, max = 10.0,),
                inputactivepowerlimit = 10.0,
                outputactivepowerlimit = 10.0,
                efficiency = (in = 0.90, out = 0.80),
                );
sys5b = PowerSystem(nodes5, generators5, loads5_DA, branches5, [battery],  100.0)
;

m = Model()

PTDF, = PowerSystems.buildptdf(sys.branches, sys.buses)
RHS = BLAS.gemm('N','N', PTDF, timeseries_netinjection)


text:
function clearswi_main(args)
    settings = getargs(args)
    if isnothing(settings) return end

    writedir = settings["output"]
    filename = "clearswi"
    if occursin(r"\.nii$", writedir)
        filename = basename(writedir)
        writedir = dirname(writedir)
    end

    mkpath(writedir)

    mag = readmag(settings["magnitude"]; mmap=!settings["no-mmap"])
    phase = readphase(settings["phase"]; mmap=!settings["no-mmap"], rescale=!settings["no-phase-rescale"])
    hdr = header(mag)
    neco = size(mag, 4)

    ## Echoes for unwrapping
    echoes = try
        getechoes(settings, neco)
    catch y
        if isa(y, BoundsError)
            error("echoes=$(join(settings["unwrap-echoes"], " ")): specified echo out of range! Number of echoes is $neco")
        else
            error("echoes=$(join(settings["unwrap-echoes"], " ")) wrongly formatted!")
        end
    end
    settings["verbose"] && println("Echoes are $echoes")

    TEs = getTEs(settings, neco, echoes)
    settings["verbose"] && println("TEs are $TEs")

    ## Error messages
    if 1 < length(echoes) && length(echoes) != length(TEs)
        error("Number of chosen echoes is $(length(echoes)) ($neco in .nii data), but $(length(TEs)) TEs were specified!")
    end

    echoes = getechoes(settings, neco)
    if echoes != 1:neco
        phase = phase[:,:,:,echoes]
        mag = mag[:,:,:,echoes]
        settings["verbose"] && println("Selecting echoes $echoes")
    end

    data = Data(mag, phase, hdr, TEs)
    mag_combine =   if settings["mag-combine"][1] == "SNR"
                        :SNR
                    elseif settings["mag-combine"][1] == "average"
                        :average
                    elseif settings["mag-combine"][1] == "echo"
                        :echo => parse(Int, last(settings["mag-combine"]))
                    elseif settings["mag-combine"][1] == "SE"
                        :SE => parse(Float32, last(settings["mag-combine"]))
                    else
                        error("The setting for mag-combine is not valid: $(settings["mag-combine"])")
                    end
    mag_sens =  if settings["mag-sensitivity-correction"] == "on"
                    nothing
                elseif settings["mag-sensitivity-correction"] == "off"
                    [1]
                elseif isfile(settings["mag-sensitivity-correction"])
                    settings["mag-sensitivity-correction"]
                else
                    error("The setting for mag-sensitivity-correction is not valid: $(settings["mag-sensitivity-correction"])")
                end
    mag_softplus =  if settings["mag-softplus-scaling"] == "on"
                        true
                    elseif settings["mag-softplus-scaling"] == "off"
                        false
                    else
                        error("The setting for mag-softplus-scaling is not valid: $(settings["mag-softplus-scaling"])")
                    end
    phase_unwrap = Symbol(settings["unwrapping-algorithm"])
    phase_hp_σ = eval(Meta.parse(join(settings["filter-size"], " ")))
    phase_scaling_type = Symbol(settings["phase-scaling-type"])
    phase_scaling_strength = try parse(Int, settings["phase-scaling-strength"]) catch; parse(Float32, settings["phase-scaling-strength"]) end
    writesteps = settings["writesteps"]

    options = Options(;mag_combine, mag_sens, mag_softplus, phase_unwrap, phase_hp_σ, phase_scaling_type, phase_scaling_strength, writesteps)

    swi = calculateSWI(data, options)
    mip = createIntensityProjection(swi, minimum)

    savenii(swi, filename, writedir, hdr)
    savenii(mip, "mip", writedir, hdr)

    return 0
end


text:
using GLMakie, Makie, Random, Colors, LinearAlgebra
using Makie.GeometryBasics: Cylinder, Pyramid
import GeometryBasics
Random.seed!(3)

cyl = Cylinder(Point{3, Float64}(1,2,3), Point{3, Float64}(2,3,4), 1.0)
pyr = Pyramid(Point3f0(0), 1f0, 1f0)
rectmesh = Rect3(Point3f(-0.5), Vec3f(1))
rectthin = Rect3(Point3f(-1), Vec3f(2,2,0.25))
sphere = Sphere(Point3f(-0.5), 1)
Cone(; quality = 10) = merge([
    Makie._circle(Point3f(0), 0.5f0, Vec3f(0,0,-1), quality),
    Makie._mantle(Point3f(0), Point3f(0,0,1), 0.5f0, 0f0, quality)])
cone = Cone()
rectMesh = GeometryBasics.mesh(rectmesh)
rectThin = GeometryBasics.mesh(rectthin)
cyL = GeometryBasics.mesh(cyl)

cmap = to_colormap(:Spectral_11, length(rectMesh.position))
colors1 = [cmap[i] for (i,v) in enumerate(rectMesh.position)]
colors2 = [RGBA(rand(4)...) for v in rectThin.position]
colors3 = [norm(v) for v in cyL.position]
markers = [sphere, rectmesh, cyl, pyr, cone]

set_theme!(theme_dark())
fig = Figure(resolution = (1200,800))
axs = [Axis3(fig[i,j]; aspect = :data, perspectiveness = 0.5)
    for j in 1:3, i in 1:2]
mesh!(axs[1], sphere, color = :white)
mesh!(axs[2], rectmesh, color = colors1)
mesh!(axs[3], pyr; color = (:dodgerblue, 0.85))
wireframe!(axs[3], pyr; color = :grey90)
mesh!(axs[4], cyl; color = colors3,
    colormap = :diverging_tritanopic_cwr_75_98_c20_n256)
mesh!(axs[5], cone; transparency = true)
wireframe!(axs[5], cone; color = :grey90, linewidth = 0.5)
mesh!(axs[6], rectthin; color = colors2, shading = false)
[meshscatter!(axs[6], Point3f(1.5rand(3) .- 0.5); marker = markers[i],
    markersize = 0.25) for i in 1:5]
save("meshes.png", fig)
display(fig)


text:
using ExactOptimalTransport
using SafeTestsets

using Test

@testset "ExactOptimalTransport" begin
    @safetestset "Utilities" begin
        include("utils.jl")
    end

    @safetestset "Exact OT" begin
        include("exact.jl")
    end

    @safetestset "Wasserstein distance" begin
        include("wasserstein.jl")
    end

    @safetestset "Bures distance" begin
        include("bures.jl")
    end
end


text:
module SessionActions

import ..Pluto: ServerSession, Notebook, emptynotebook, tamepath, move_notebook!, update_save_run!, putnotebookupdates!, putplutoupdates!, load_notebook, get_pl_env, clientupdate_notebook_list, WorkspaceManager, @asynclog

struct NotebookIsRunningException <: Exception
    notebook::Notebook
end

function open_url(session::ServerSession, url::AbstractString; kwargs...)
    path = download(url, emptynotebook().path)
    open(session, path; kwargs...)
end

function open(session::ServerSession, path::AbstractString; run_async=true)
    for nb in values(session.notebooks)
        if realpath(nb.path) == realpath(tamepath(path))
            throw(NotebookIsRunningException(nb))
        end
    end

    nb = load_notebook(tamepath(path))
    session.notebooks[nb.notebook_id] = nb
    if get_pl_env("PLUTO_RUN_NOTEBOOK_ON_LOAD") == "true"
        update_save_run!(session, nb, nb.cells; run_async=run_async, prerender_text=true)
        # TODO: send message when initial run completed
    end

    if run_async
        @asynclog putplutoupdates!(session, clientupdate_notebook_list(session.notebooks))
    else
        putplutoupdates!(session, clientupdate_notebook_list(session.notebooks))
    end

    nb
end

function new(session::ServerSession; run_async=true)
    nb = emptynotebook()
    update_save_run!(session, nb, nb.cells; run_async=run_async, prerender_text=true)
    session.notebooks[nb.notebook_id] = nb

    if run_async
        @asynclog putplutoupdates!(session, clientupdate_notebook_list(session.notebooks))
    else
        putplutoupdates!(session, clientupdate_notebook_list(session.notebooks))
    end

    nb
end

function shutdown(session::ServerSession, notebook::Notebook; keep_in_session=false)
    if !keep_in_session
        listeners = putnotebookupdates!(session, notebook) # TODO: shutdown message
        delete!(session.notebooks, notebook.notebook_id)
        putplutoupdates!(session, clientupdate_notebook_list(session.notebooks))
        for client in listeners
            @async close(client.stream)
        end
    end
    success = WorkspaceManager.unmake_workspace(notebook)
end

function move(session::ServerSession, notebook::Notebook, newpath::AbstractString)
    result = try
        if isfile(newpath)
            (success = false, reason = "File exists already - you need to delete the old file manually.")
        else
            move_notebook!(notebook, newpath)
            putplutoupdates!(session, clientupdate_notebook_list(session.notebooks))
            WorkspaceManager.cd_workspace(notebook, newpath)
            (success = true, reason = "")
        end
    catch ex
        showerror(stderr, stacktrace(catch_backtrace()))
        (success = false, reason = sprint(showerror, ex))
    end
    result
end
end


text:
@doc raw"""
    IdealGlmMhdEquations1D

The ideal compressible GLM-MHD equations in one space dimension.
!!! note
    There is no divergence cleaning variable `psi` because the divergence-free constraint
    is satisfied trivially in one spatial dimension.
"""
struct IdealGlmMhdEquations1D{RealT<:Real} <: AbstractIdealGlmMhdEquations{1, 8}
  gamma::RealT
end

have_nonconservative_terms(::IdealGlmMhdEquations1D) = Val(false)
varnames(::typeof(cons2cons), ::IdealGlmMhdEquations1D) = ("rho", "rho_v1", "rho_v2", "rho_v3", "rho_e", "B1", "B2", "B3")
varnames(::typeof(cons2prim), ::IdealGlmMhdEquations1D) = ("rho", "v1", "v2", "v3", "p", "B1", "B2", "B3")
default_analysis_integrals(::IdealGlmMhdEquations1D)  = (entropy_timederivative, Val(:l2_divb), Val(:linf_divb))

"""
    initial_condition_constant(x, t, equations::IdealGlmMhdEquations1D)

A constant initial condition to test free-stream preservation.
"""
function initial_condition_constant(x, t, equations::IdealGlmMhdEquations1D)
  rho = 1.0
  rho_v1 = 0.1
  rho_v2 = -0.2
  rho_v3 = -0.5
  rho_e = 50.0
  B1 = 3.0
  B2 = -1.2
  B3 = 0.5
  return SVector(rho, rho_v1, rho_v2, rho_v3, rho_e, B1, B2, B3)
end

"""
    initial_condition_convergence_test(x, t, equations::IdealGlmMhdEquations1D)

An Alfvén wave as smooth initial condition used for convergence tests.
"""
function initial_condition_convergence_test(x, t, equations::IdealGlmMhdEquations1D)
  # smooth Alfvén wave test from Derigs et al. FLASH (2016)
  # domain must be set to [0, 1], γ = 5/3
  rho = 1.0
  v1 = 0.0
  # TODO: sincospi
  si, co = sincos(2 * pi * x[1])
  v2 = 0.1 * si
  v3 = 0.1 * co
  p = 0.1
  B1 = 1.0
  B2 = v2
  B3 = v3
  return prim2cons(SVector(rho, v1, v2, v3, p, B1, B2, B3), equations)
end

"""
    initial_condition_briowu_shock_tube(x, t, equations::IdealGlmMhdEquations1D)

Compound shock tube test case for one dimensional ideal MHD equations. It is bascially an
MHD extension of the Sod shock tube. Taken from Section V of the article
- Brio and Wu (1988)
  An Upwind Differencing Scheme for the Equations of Ideal Magnetohydrodynamics
  [DOI: 10.1016/0021-9991(88)90120-9](https://doi.org/10.1016/0021-9991(88)90120-9)
"""
function initial_condition_briowu_shock_tube(x, t, equations::IdealGlmMhdEquations1D)
  # domain must be set to [0, 1], γ = 2, final time = 0.12
  rho = x[1] < 0.5 ? 1.0 : 0.125
  v1 = 0.0
  v2 = 0.0
  v3 = 0.0
  p = x[1] < 0.5 ? 1.0 : 0.1
  B1 = 0.75
  B2 = x[1] < 0.5 ? 1.0 : -1.0
  B3 = 0.0
  return prim2cons(SVector(rho, v1, v2, v3, p, B1, B2, B3), equations)
end

"""
    boundary_condition_briowu_shock_tube(u_inner, orientation, direction, x, t,
                                         surface_flux_function,
                                         equations::IdealGlmMhdEquations1D)

Boundary conditions used for the Brio and Wu shock tube in combination with
[`initial_condition_briowu_shock_tube`](@ref).
"""
function boundary_condition_briowu_shock_tube(u_inner, orientation, direction, x, t,
                                              surface_flux_function,
                                              equations::IdealGlmMhdEquations1D)
  u_boundary = initial_condition_briowu_shock_tube(x, t, equations)
  # Calculate boundary flux
  if direction == 2 # u_inner is "left" of boundary, u_boundary is "right" of boundary
    flux = surface_flux_function(u_inner, u_boundary, orientation, equations)
  else # u_boundary is "left" of boundary, u_inner is "right" of boundary
    flux = surface_flux_function(u_boundary, u_inner, orientation, equations)
  end

  return flux
end

"""
    initial_condition_torrilhon_shock_tube(x, t, equations::IdealGlmMhdEquations1D)

Torrilhon's shock tube test case for one dimensional ideal MHD equations.
- Torrilhon (2003)
  Uniqueness conditions for Riemann problems of ideal magnetohydrodynamics
  [DOI: 10.1017/S0022377803002186](https://doi.org/10.1017/S0022377803002186)
"""
function initial_condition_torrilhon_shock_tube(x, t, equations::IdealGlmMhdEquations1D)
  # domain must be set to [-1, 1.5], γ = 5/3, final time = 0.4
  rho = x[1] <= 0 ? 3.0 : 1.0
  v1 = 0.0
  v2 = 0.0
  v3 = 0.0
  p = x[1] <= 0 ? 3.0 : 1.0
  B1 = 1.5
  B2 = x[1] <= 0 ? 1.0 : cos(1.5)
  B3 = x[1] <= 0 ? 0.0 : sin(1.5)
  return prim2cons(SVector(rho, v1, v2, v3, p, B1, B2, B3), equations)
end

"""
    boundary_condition_torrilhon_shock_tube(u_inner, orientation, direction, x, t,
                                            surface_flux_function,
                                            equations::IdealGlmMhdEquations1D)

Boundary conditions used for the Torrilhon shock tube in combination with
[`initial_condition_torrilhon_shock_tube`](@ref).
"""
function boundary_condition_torrilhon_shock_tube(u_inner, orientation, direction, x, t,
                                                 surface_flux_function,
                                                 equations::IdealGlmMhdEquations1D)
  u_boundary = initial_condition_torrilhon_shock_tube(x, t, equations)
  # Calculate boundary flux
  if direction == 2 # u_inner is "left" of boundary, u_boundary is "right" of boundary
    flux = surface_flux_function(u_inner, u_boundary, orientation, equations)
  else # u_boundary is "left" of boundary, u_inner is "right" of boundary
    flux = surface_flux_function(u_boundary, u_inner, orientation, equations)
  end

  return flux
end

"""
    initial_condition_ryujones_shock_tube(x, t, equations::IdealGlmMhdEquations1D)

Ryu and Jones shock tube test case for one dimensional ideal MHD equations. Contains
fast shocks, slow shocks, and rational discontinuities that propagate on either side
of the contact discontinuity. Exercises the scheme to capture all 7 types of waves
present in the one dimensional MHD equations. It is the second test from Section 4 of
- Ryu and Jones (1995)
  Numerical Magnetohydrodynamics in Astrophysics: Algorithm and Tests
  for One-Dimensional Flow
  [DOI: 10.1086/175437](https://doi.org/10.1086/175437)
!!! note
    This paper has a typo in the initial conditions. Their variable `E` should be `p`.
"""
function initial_condition_ryujones_shock_tube(x, t, equations::IdealGlmMhdEquations1D)
  # domain must be set to [0, 1], γ = 5/3, final time = 0.2
  rho = x[1] <= 0.5 ? 1.08 : 1.0
  v1 = x[1] <= 0.5 ? 1.2 : 0.0
  v2 = x[1] <= 0.5 ? 0.01 : 0.0
  v3 = x[1] <= 0.5 ? 0.5 : 0.0
  p = x[1] <= 0.5 ? 0.95 : 1.0
  inv_sqrt4pi = 1.0 / sqrt(4 * pi)
  B1 = 2 * inv_sqrt4pi
  B2 = x[1] <= 0.5 ? 3.6 * inv_sqrt4pi : 4.0 * inv_sqrt4pi
  B3 = B1

  return prim2cons(SVector(rho, v1, v2, v3, p, B1, B2, B3), equations)
end

"""
    boundary_condition_ryujones_shock_tube(u_inner, orientation, direction, x, t,
                                            surface_flux_function,
                                            equations::IdealGlmMhdEquations1D)

Boundary conditions used for the Ryu and Jones shock tube in combination with
[`initial_condition_ryujones_shock_tube`](@ref).
"""
function boundary_condition_ryujones_shock_tube(u_inner, orientation, direction, x, t,
                                                surface_flux_function,
                                                equations::IdealGlmMhdEquations1D)
  u_boundary = initial_condition_ryujones_shock_tube(x, t, equations)
  # Calculate boundary flux
  if direction == 2 # u_inner is "left" of boundary, u_boundary is "right" of boundary
    flux = surface_flux_function(u_inner, u_boundary, orientation, equations)
  else # u_boundary is "left" of boundary, u_inner is "right" of boundary
    flux = surface_flux_function(u_boundary, u_inner, orientation, equations)
  end

  return flux
end

"""
    initial_condition_shu_osher_shock_tube(x, t, equations::IdealGlmMhdEquations1D)

Extended version of the test of Shu and Osher for one dimensional ideal MHD equations.
Taken from Section 4.1 of
- Derigs et al. (2016)
  A Novel High-Order, Entropy Stable, 3D AMR MHD Solver withGuaranteed Positive Pressure
  [DOI: 10.1016/j.jcp.2016.04.048](https://doi.org/10.1016/j.jcp.2016.04.048)
"""
function initial_condition_shu_osher_shock_tube(x, t, equations::IdealGlmMhdEquations1D)
  # domain must be set to [-5, 5], γ = 5/3, final time = 0.7
  # initial shock location is taken to be at x = -4
  x_0 = -4.0
  rho = x[1] <= x_0 ? 3.5 : 1.0 + 0.2 * sin(5.0 * x[1])
  v1 = x[1] <= x_0 ? 5.8846 : 0.0
  v2 = x[1] <= x_0 ? 1.1198 : 0.0
  v3 = 0.0
  p = x[1] <= x_0 ? 42.0267 : 1.0
  B1 = 1.0
  B2 = x[1] <= x_0 ? 3.6359 : 1.0
  B3 = 0.0

  return prim2cons(SVector(rho, v1, v2, v3, p, B1, B2, B3), equations)
end

"""
    boundary_condition_shu_osher_shock_tube(u_inner, orientation, direction, x, t,
                                            surface_flux_function,
                                            equations::IdealGlmMhdEquations1D)

Boundary conditions used for the MHD extension of the Shu-Osher test case
in combination with [`initial_condition_shu_osher_shock_tube`](@ref).
"""
function boundary_condition_shu_osher_shock_tube(u_inner, orientation, direction, x, t,
                                                surface_flux_function,
                                                equations::IdealGlmMhdEquations1D)
  u_boundary = initial_condition_shu_osher_shock_tube(x, t, equations)
  # Calculate boundary flux
  if direction == 2 # u_inner is "left" of boundary, u_boundary is "right" of boundary
    flux = surface_flux_function(u_inner, u_boundary, orientation, equations)
  else # u_boundary is "left" of boundary, u_inner is "right" of boundary
    flux = surface_flux_function(u_boundary, u_inner, orientation, equations)
  end

  return flux
end

"""
    initial_condition_shu_osher_shock_tube_flipped(x, t, equations::IdealGlmMhdEquations1D)

Extended version of the test of Shu and Osher for one dimensional ideal MHD equations
but shock propogates from right to left.
!!! note
    This is useful to exercise some of the components of the HLL flux.
"""
function initial_condition_shu_osher_shock_tube_flipped(x, t, equations::IdealGlmMhdEquations1D)
  # domain must be set to [-5, 5], γ = 5/3, final time = 0.7
  # initial shock location is taken to be at x = 4
  x_0 = 4.0
  rho = x[1] <= x_0 ? 1.0 + 0.2 * sin(5.0 * x[1]) : 3.5
  v1 = x[1] <= x_0 ? 0.0 : -5.8846
  v2 = x[1] <= x_0 ? 0.0 : -1.1198
  v3 = 0.0
  p = x[1] <= x_0 ? 1.0 : 42.0267
  B1 = 1.0
  B2 = x[1] <= x_0 ? 1.0 : 3.6359
  B3 = 0.0

  return prim2cons(SVector(rho, v1, v2, v3, p, B1, B2, B3), equations)
end

"""
    boundary_condition_shu_osher_shock_tube_flipped(u_inner, orientation, direction, x, t,
                                            surface_flux_function,
                                            equations::IdealGlmMhdEquations1D)

Boundary conditions used for the MHD extension of the Shu-Osher test case that travels
from right to left. Use with [`initial_condition_shu_osher_shock_tube_flipped`](@ref).
"""
function boundary_condition_shu_osher_shock_tube_flipped(u_inner, orientation, direction, x, t,
                                                surface_flux_function,
                                                equations::IdealGlmMhdEquations1D)
  u_boundary = initial_condition_shu_osher_shock_tube_flipped(x, t, equations)
  # Calculate boundary flux
  if direction == 2 # u_inner is "left" of boundary, u_boundary is "right" of boundary
    flux = surface_flux_function(u_inner, u_boundary, orientation, equations)
  else # u_boundary is "left" of boundary, u_inner is "right" of boundary
    flux = surface_flux_function(u_boundary, u_inner, orientation, equations)
  end

  return flux
end

# Calculate 1D flux in for a single point
@inline function flux(u, orientation::Integer, equations::IdealGlmMhdEquations1D)
  rho, rho_v1, rho_v2, rho_v3, rho_e, B1, B2, B3 = u
  v1 = rho_v1/rho
  v2 = rho_v2/rho
  v3 = rho_v3/rho
  kin_en = 0.5 * rho * (v1^2 + v2^2 + v3^2)
  mag_en = 0.5 * (B1^2 + B2^2 + B3^2)
  p = (equations.gamma - 1) * (rho_e - kin_en - mag_en)

  # Ignore orientation since it is always "1" in 1D
  f1 = rho_v1
  f2 = rho_v1*v1 + p + mag_en - B1^2
  f3 = rho_v1*v2 - B1*B2
  f4 = rho_v1*v3 - B1*B3
  f5 = (kin_en + equations.gamma*p/(equations.gamma - 1) + 2*mag_en)*v1 - B1*(v1*B1 + v2*B2 + v3*B3)
  f6 = 0.0
  f7 = v1*B2 - v2*B1
  f8 = v1*B3 - v3*B1

  return SVector(f1, f2, f3, f4, f5, f6, f7, f8)
end

"""
    flux_derigs_etal(u_ll, u_rr, orientation, equations::IdealGlmMhdEquations1D)

Entropy conserving two-point flux by
- Derigs et al. (2018)
  Ideal GLM-MHD: About the entropy consistent nine-wave magnetic field
  divergence diminishing ideal magnetohydrodynamics equations
  [DOI: 10.1016/j.jcp.2018.03.002](https://doi.org/10.1016/j.jcp.2018.03.002)
"""
function flux_derigs_etal(u_ll, u_rr, orientation::Integer, equations::IdealGlmMhdEquations1D)
  # Unpack left and right states to get velocities, pressure, and inverse temperature (called beta)
  rho_ll, rho_v1_ll, rho_v2_ll, rho_v3_ll, rho_e_ll, B1_ll, B2_ll, B3_ll = u_ll
  rho_rr, rho_v1_rr, rho_v2_rr, rho_v3_rr, rho_e_rr, B1_rr, B2_rr, B3_rr = u_rr

  v1_ll = rho_v1_ll/rho_ll
  v2_ll = rho_v2_ll/rho_ll
  v3_ll = rho_v3_ll/rho_ll
  v1_rr = rho_v1_rr/rho_rr
  v2_rr = rho_v2_rr/rho_rr
  v3_rr = rho_v3_rr/rho_rr
  vel_norm_ll = v1_ll^2 + v2_ll^2 + v3_ll^2
  vel_norm_rr = v1_rr^2 + v2_rr^2 + v3_rr^2
  mag_norm_ll = B1_ll^2 + B2_ll^2 + B3_ll^2
  mag_norm_rr = B1_rr^2 + B2_rr^2 + B3_rr^2
  p_ll = (equations.gamma - 1)*(rho_e_ll - 0.5*rho_ll*vel_norm_ll - 0.5*mag_norm_ll)
  p_rr = (equations.gamma - 1)*(rho_e_rr - 0.5*rho_rr*vel_norm_rr - 0.5*mag_norm_rr)
  beta_ll = 0.5*rho_ll/p_ll
  beta_rr = 0.5*rho_rr/p_rr
  # for convenience store v⋅B
  vel_dot_mag_ll = v1_ll*B1_ll + v2_ll*B2_ll + v3_ll*B3_ll
  vel_dot_mag_rr = v1_rr*B1_rr + v2_rr*B2_rr + v3_rr*B3_rr

  # Compute the necessary mean values needed for either direction
  rho_avg  = 0.5*(rho_ll+rho_rr)
  rho_mean = ln_mean(rho_ll,rho_rr)
  beta_mean = ln_mean(beta_ll,beta_rr)
  beta_avg = 0.5*(beta_ll+beta_rr)
  v1_avg = 0.5*(v1_ll+v1_rr)
  v2_avg = 0.5*(v2_ll+v2_rr)
  v3_avg = 0.5*(v3_ll+v3_rr)
  p_mean = 0.5*rho_avg/beta_avg
  B1_avg = 0.5*(B1_ll+B1_rr)
  B2_avg = 0.5*(B2_ll+B2_rr)
  B3_avg = 0.5*(B3_ll+B3_rr)
  vel_norm_avg = 0.5*(vel_norm_ll+vel_norm_rr)
  mag_norm_avg = 0.5*(mag_norm_ll+mag_norm_rr)
  vel_dot_mag_avg = 0.5*(vel_dot_mag_ll+vel_dot_mag_rr)

  # Ignore orientation since it is always "1" in 1D
  f1 = rho_mean*v1_avg
  f2 = f1*v1_avg + p_mean + 0.5*mag_norm_avg - B1_avg*B1_avg
  f3 = f1*v2_avg - B1_avg*B2_avg
  f4 = f1*v3_avg - B1_avg*B3_avg
  f6 = 0.0
  f7 = v1_avg*B2_avg - v2_avg*B1_avg
  f8 = v1_avg*B3_avg - v3_avg*B1_avg
  # total energy flux is complicated and involves the previous eight components
  v1_mag_avg = 0.5*(v1_ll*mag_norm_ll + v1_rr*mag_norm_rr)
  f5 = (f1*0.5*(1/(equations.gamma-1)/beta_mean - vel_norm_avg) + f2*v1_avg + f3*v2_avg +
        f4*v3_avg + f6*B1_avg + f7*B2_avg + f8*B3_avg - 0.5*v1_mag_avg +
        B1_avg*vel_dot_mag_avg)

  return SVector(f1, f2, f3, f4, f5, f6, f7, f8)
end

# Calculate maximum wave speed for local Lax-Friedrichs-type dissipation
@inline function max_abs_speed_naive(u_ll, u_rr, orientation::Integer, equations::IdealGlmMhdEquations1D)
  rho_ll, rho_v1_ll, rho_v2_ll, rho_v3_ll, rho_e_ll, B1_ll, B2_ll, B3_ll = u_ll
  rho_rr, rho_v1_rr, rho_v2_rr, rho_v3_rr, rho_e_rr, B1_rr, B2_rr, B3_rr = u_rr

  # Calculate velocities and fast magnetoacoustic wave speeds
  # left
  v1_ll = rho_v1_ll / rho_ll
  v2_ll = rho_v2_ll / rho_ll
  v3_ll = rho_v3_ll / rho_ll
  v_mag_ll = sqrt(v1_ll^2 + v2_ll^2 + v3_ll^2)
  cf_ll = calc_fast_wavespeed(u_ll, orientation, equations)
  # right
  v1_rr = rho_v1_rr / rho_rr
  v2_rr = rho_v2_rr / rho_rr
  v3_rr = rho_v3_rr / rho_rr
  v_mag_rr = sqrt(v1_rr^2 + v2_rr^2 + v3_rr^2)
  cf_rr = calc_fast_wavespeed(u_rr, orientation, equations)

  λ_max = max(v_mag_ll, v_mag_rr) + max(cf_ll, cf_rr)
end

"""
    min_max_speed_naive(u_ll, u_rr, orientation, equations::IdealGlmMhdEquations1D)

Calculate minimum and maximum wave speeds for HLL-type fluxes as in
- Li (2005)
  An HLLC Riemann solver for magneto-hydrodynamics
  [DOI: 10.1016/j.jcp.2004.08.020](https://doi.org/10.1016/j.jcp.2004.08.020)
"""
@inline function min_max_speed_naive(u_ll, u_rr, orientation::Integer, equations::IdealGlmMhdEquations1D)
  rho_ll, rho_v1_ll, rho_v2_ll, rho_v3_ll, rho_e_ll, B1_ll, B2_ll, B3_ll = u_ll
  rho_rr, rho_v1_rr, rho_v2_rr, rho_v3_rr, rho_e_rr, B1_rr, B2_rr, B3_rr = u_rr

  # Calculate primitive variables and speed of sound
  v1_ll = rho_v1_ll/rho_ll
  v2_ll = rho_v2_ll/rho_ll
  v3_ll = rho_v3_ll/rho_ll
  vel_norm_ll = v1_ll^2 + v2_ll^2 + v3_ll^2
  mag_norm_ll = B1_ll^2 + B2_ll^2 + B3_ll^2
  p_ll = (equations.gamma - 1)*(rho_e_ll - 0.5*rho_ll*vel_norm_ll - 0.5*mag_norm_ll)

  v1_rr = rho_v1_rr/rho_rr
  v2_rr = rho_v2_rr/rho_rr
  v3_rr = rho_v3_rr/rho_rr
  vel_norm_rr = v1_rr^2 + v2_rr^2 + v3_rr^2
  mag_norm_rr = B1_rr^2 + B2_rr^2 + B3_rr^2
  p_rr = (equations.gamma - 1)*(rho_e_rr - 0.5*rho_rr*vel_norm_rr - 0.5*mag_norm_rr)

  # Approximate the left-most and right-most eigenvalues in the Riemann fan
  c_f_ll = calc_fast_wavespeed(u_ll, orientation, equations)
  c_f_rr = calc_fast_wavespeed(u_rr, orientation, equations)
  vel_roe, c_f_roe = calc_fast_wavespeed_roe(u_ll, u_rr, orientation, equations)
  λ_min = min(v1_ll - c_f_ll, vel_roe - c_f_roe)
  λ_max = max(v1_rr + c_f_rr, vel_roe + c_f_roe)

  return λ_min, λ_max
end

@inline function max_abs_speeds(u, equations::IdealGlmMhdEquations1D)
  rho, rho_v1, _ = u
  v1 = rho_v1 / rho
  cf_x_direction = calc_fast_wavespeed(u, 1, equations)

  return abs(v1) + cf_x_direction
end

# Convert conservative variables to primitive
function cons2prim(u, equations::IdealGlmMhdEquations1D)
  rho, rho_v1, rho_v2, rho_v3, rho_e, B1, B2, B3 = u

  v1 = rho_v1 / rho
  v2 = rho_v2 / rho
  v3 = rho_v3 / rho
  p = (equations.gamma - 1) * (rho_e - 0.5*rho*(v1^2 + v2^2 + v3^2) - 0.5*(B1^2 + B2^2 + B3^2))

  return SVector(rho, v1, v2, v3, p, B1, B2, B3)
end

# Convert conservative variables to entropy
@inline function cons2entropy(u, equations::IdealGlmMhdEquations1D)
  rho, rho_v1, rho_v2, rho_v3, rho_e, B1, B2, B3 = u

  v1 = rho_v1 / rho
  v2 = rho_v2 / rho
  v3 = rho_v3 / rho
  v_square = v1^2 + v2^2 + v3^2
  p = (equations.gamma - 1) * (rho_e - 0.5*rho*v_square - 0.5*(B1^2 + B2^2 + B3^2))
  s = log(p) - equations.gamma*log(rho)
  rho_p = rho / p

  w1 = (equations.gamma - s) / (equations.gamma-1) - 0.5 * rho_p * v_square
  w2 = rho_p * v1
  w3 = rho_p * v2
  w4 = rho_p * v3
  w5 = -rho_p
  w6 = rho_p * B1
  w7 = rho_p * B2
  w8 = rho_p * B3

  return SVector(w1, w2, w3, w4, w5, w6, w7, w8)
end

# Convert primitive to conservative variables
@inline function prim2cons(prim, equations::IdealGlmMhdEquations1D)
  rho, v1, v2, v3, p, B1, B2, B3 = prim

  rho_v1 = rho * v1
  rho_v2 = rho * v2
  rho_v3 = rho * v3
  rho_e = p/(equations.gamma-1) + 0.5 * (rho_v1*v1 + rho_v2*v2 + rho_v3*v3) +
                                 0.5 * (B1^2 + B2^2 + B3^2)

  return SVector(rho, rho_v1, rho_v2, rho_v3, rho_e, B1, B2, B3)
end

@inline function density(u, equations::IdealGlmMhdEquations1D)
  rho, rho_v1, rho_v2, rho_v3, rho_e, B1, B2, B3 = u
  return rho
end

@inline function pressure(u, equations::IdealGlmMhdEquations1D)
  rho, rho_v1, rho_v2, rho_v3, rho_e, B1, B2, B3 = u
  p = (equations.gamma - 1)*(rho_e - 0.5 * (rho_v1^2 + rho_v2^2 + rho_v3^2) / rho
                                   - 0.5 * (B1^2 + B2^2 + B3^2))
  return p
end

@inline function density_pressure(u, equations::IdealGlmMhdEquations1D)
  rho, rho_v1, rho_v2, rho_v3, rho_e, B1, B2, B3 = u
  p = (equations.gamma - 1)*(rho_e - 0.5 * (rho_v1^2 + rho_v2^2 + rho_v3^2) / rho
                                   - 0.5 * (B1^2 + B2^2 + B3^2))
  return rho * p
end

# Compute the fastest wave speed for ideal MHD equations: c_f, the fast magnetoacoustic eigenvalue
@inline function calc_fast_wavespeed(cons, direction, equations::IdealGlmMhdEquations1D)
  rho, rho_v1, rho_v2, rho_v3, rho_e, B1, B2, B3 = cons
  v1 = rho_v1 / rho
  v2 = rho_v2 / rho
  v3 = rho_v3 / rho
  v_mag = sqrt(v1^2 + v2^2 + v3^2)
  p = (equations.gamma - 1)*(rho_e - 0.5*rho*v_mag^2 - 0.5*(B1^2 + B2^2 + B3^2))
  a_square = equations.gamma * p / rho
  sqrt_rho = sqrt(rho)
  b1 = B1 / sqrt_rho
  b2 = B2 / sqrt_rho
  b3 = B3 / sqrt_rho
  b_square = b1^2 + b2^2 + b3^2

  c_f = sqrt(0.5*(a_square + b_square) + 0.5*sqrt((a_square + b_square)^2 - 4.0*a_square*b1^2))
  return c_f
end

"""
    calc_fast_wavespeed_roe(u_ll, u_rr, direction, equations::IdealGlmMhdEquations1D)

Compute the fast magnetoacoustic wave speed using Roe averages
as given by
- Cargo and Gallice (1997)
  Roe Matrices for Ideal MHD and Systematic Construction
  of Roe Matrices for Systems of Conservation Laws
  [DOI: 10.1006/jcph.1997.5773](https://doi.org/10.1006/jcph.1997.5773)
"""
@inline function calc_fast_wavespeed_roe(u_ll, u_rr, direction, equations::IdealGlmMhdEquations1D)
  rho_ll, rho_v1_ll, rho_v2_ll, rho_v3_ll, rho_e_ll, B1_ll, B2_ll, B3_ll = u_ll
  rho_rr, rho_v1_rr, rho_v2_rr, rho_v3_rr, rho_e_rr, B1_rr, B2_rr, B3_rr = u_rr

  # Calculate primitive variables
  v1_ll = rho_v1_ll/rho_ll
  v2_ll = rho_v2_ll/rho_ll
  v3_ll = rho_v3_ll/rho_ll
  vel_norm_ll = v1_ll^2 + v2_ll^2 + v3_ll^2
  mag_norm_ll = B1_ll^2 + B2_ll^2 + B3_ll^2
  p_ll = (equations.gamma - 1)*(rho_e_ll - 0.5*rho_ll*vel_norm_ll - 0.5*mag_norm_ll)

  v1_rr = rho_v1_rr/rho_rr
  v2_rr = rho_v2_rr/rho_rr
  v3_rr = rho_v3_rr/rho_rr
  vel_norm_rr = v1_rr^2 + v2_rr^2 + v3_rr^2
  mag_norm_rr = B1_rr^2 + B2_rr^2 + B3_rr^2
  p_rr = (equations.gamma - 1)*(rho_e_rr - 0.5*rho_rr*vel_norm_rr - 0.5*mag_norm_rr)

  # compute total pressure which is thermal + magnetic pressures
  p_total_ll = p_ll + 0.5*mag_norm_ll
  p_total_rr = p_rr + 0.5*mag_norm_rr

  # compute the Roe density averages
  sqrt_rho_ll = sqrt(rho_ll)
  sqrt_rho_rr = sqrt(rho_rr)
  inv_sqrt_rho_add  = 1.0 / (sqrt_rho_ll + sqrt_rho_rr)
  inv_sqrt_rho_prod = 1.0 / (sqrt_rho_ll * sqrt_rho_rr)
  rho_ll_roe =  sqrt_rho_ll * inv_sqrt_rho_add
  rho_rr_roe =  sqrt_rho_rr * inv_sqrt_rho_add
  # Roe averages
  # velocities and magnetic fields
  v1_roe = v1_ll * rho_ll_roe + v1_rr * rho_rr_roe
  v2_roe = v2_ll * rho_ll_roe + v2_rr * rho_rr_roe
  v3_roe = v3_ll * rho_ll_roe + v3_rr * rho_rr_roe
  B1_roe = B1_ll * rho_ll_roe + B1_rr * rho_rr_roe
  B2_roe = B2_ll * rho_ll_roe + B2_rr * rho_rr_roe
  B3_roe = B3_ll * rho_ll_roe + B3_rr * rho_rr_roe
  # enthalpy
  H_ll  = (rho_e_ll + p_total_ll) / rho_ll
  H_rr  = (rho_e_rr + p_total_rr) / rho_rr
  H_roe = H_ll * rho_ll_roe + H_rr * rho_rr_roe
  # temporary vairable see equations (4.12) in Cargo and Gallice
  X = 0.5 * ( (B1_ll - B1_rr)^2 + (B2_ll - B2_rr)^2 + (B3_ll - B3_rr)^2 ) * inv_sqrt_rho_add^2
  # averaged components needed to compute c_f, the fast magnetoacoustic wave speed
  b_square_roe = (B1_roe^2 + B2_roe^2 + B3_roe^2) * inv_sqrt_rho_prod # scaled magnectic sum
  a_square_roe = ((2.0 - equations.gamma) * X +
                 (equations.gamma -1.0) * (H_roe - 0.5*(v1_roe^2 + v2_roe^2 + v3_roe^2) -
                                          b_square_roe)) # acoustic speed
  # finally compute the average wave speed and set the output velocity
  # Ignore orientation since it is always "1" in 1D
  c_a_roe = B1_roe^2 * inv_sqrt_rho_prod # (squared) Alfvén wave speed
  a_star_roe = sqrt( (a_square_roe + b_square_roe)^2 - 4.0 * a_square_roe * c_a_roe )
  c_f_roe = sqrt( 0.5 * (a_square_roe + b_square_roe + a_star_roe) )

  return v1_roe, c_f_roe
end

# Calculate thermodynamic entropy for a conservative state `cons`
@inline function entropy_thermodynamic(cons, equations::IdealGlmMhdEquations1D)
  # Pressure
  p = (equations.gamma - 1) * (cons[5] - 1/2 * (cons[2]^2 + cons[3]^2 + cons[4]^2) / cons[1]
                                       - 1/2 * (cons[6]^2 + cons[7]^2 + cons[8]^2))

  # Thermodynamic entropy
  s = log(p) - equations.gamma*log(cons[1])

  return s
end

# Calculate mathematical entropy for a conservative state `cons`
@inline function entropy_math(cons, equations::IdealGlmMhdEquations1D)
  S = -entropy_thermodynamic(cons, equations) * cons[1] / (equations.gamma - 1)

  return S
end

# Default entropy is the mathematical entropy
@inline entropy(cons, equations::IdealGlmMhdEquations1D) = entropy_math(cons, equations)

# Calculate total energy for a conservative state `cons`
@inline energy_total(cons, ::IdealGlmMhdEquations1D) = cons[5]

# Calculate kinetic energy for a conservative state `cons`
@inline function energy_kinetic(cons, equations::IdealGlmMhdEquations1D)
  return 0.5 * (cons[2]^2 + cons[3]^2 + cons[4]^2)/cons[1]
end

# Calculate the magnetic energy for a conservative state `cons'.
#  OBS! For non-dinmensional form of the ideal MHD magnetic pressure ≡ magnetic energy
@inline function energy_magnetic(cons, ::IdealGlmMhdEquations1D)
  return 0.5 * (cons[6]^2 + cons[7]^2 + cons[8]^2)
end

# Calculate internal energy for a conservative state `cons`
@inline function energy_internal(cons, equations::IdealGlmMhdEquations1D)
  return (energy_total(cons, equations)
          - energy_kinetic(cons, equations)
          - energy_magnetic(cons, equations))
end

# Calcluate the cross helicity (\vec{v}⋅\vec{B}) for a conservative state `cons'
@inline function cross_helicity(cons, ::IdealGlmMhdEquations1D)
  return (cons[2]*cons[6] + cons[3]*cons[7] + cons[4]*cons[8]) / cons[1]
end


text:
include("math.jl")
using KitBase

quadratureorder =2
points, triangulation = KitBase.octa_quadrature(quadratureorder)

println("-------------------")
#println(points)
println(size(points))

#nq = size(points, 1)

momentBasis =  ComputeSphericalBasisKarth(points,2, 3)
#momentBasis2 = ComputeSphericalBasisAnalytical(points)
println(size(momentBasis))
println(momentBasis)
#println(size(momentBasis2))
#println(momentBasis2-momentBasis)

#create points
points1D = linspace(-1,1,5)
points1DM = zeros(5,1)
points1DM[:,1] = points1D[:]
#println(size(points1DM))

momentBasis =  ComputeSphericalBasisKarth(points1DM,2,1)
momentBasis2 =  ComputeSphericalBasisAnalytical(points1DM)
println(size(momentBasis))
println(size(momentBasis2))
println(momentBasis-momentBasis2)


text:
using Clang

using HELICS_jll

# LIBCLANG_HEADERS are those headers to be wrapped.
const LIBHELICS_INCLUDE = [
                           normpath(joinpath(HELICS_jll.artifact_dir, "include", "helics", "shared_api_library")),
                           normpath(joinpath(HELICS_jll.artifact_dir, "include", "helics")),
                          ]

const LIBHELICS_HEADERS = String[]

for folder in LIBHELICS_INCLUDE
    for header in readdir(folder)
        if endswith(header, ".h")
            push!(LIBHELICS_HEADERS, joinpath(folder, header))
        end
    end
end

LIBHELICS_ARGS = String[]
for path in LIBHELICS_INCLUDE
    push!(LIBHELICS_ARGS, "-I")
    push!(LIBHELICS_ARGS, path)
end

wc = init(; headers = LIBHELICS_HEADERS,
            output_file = joinpath(@__DIR__, "..", "src", "lib.jl"),
            common_file = joinpath(@__DIR__, "..", "src", "common.jl"),
            clang_includes = vcat(LIBHELICS_INCLUDE..., CLANG_INCLUDE),
            clang_args = LIBHELICS_ARGS,
            header_wrapped = (root, current)->root == current,
            header_library = x->"libhelicsSharedLib",
            clang_diagnostics = true,
            )

run(wc)


text:
"""
This module deals with controlling what is displayed on the screen.

The two major features of this module are switching to the alernative screen and clearing the what is currenlty on the screen.

The [alternative screen](https://invisible-island.net/xterm/ctlseqs/ctlseqs.html#h2-The-Alternate-Screen-Buffer) is a feature of most terminals that allows an application to switch to a screen that has no scrollback, allowing it to draw whatever it wants to the terminal without interfering with the scrollback of the normal screen. This is most useful for full-screen applications, such as vim or emacs. Ansillary allows you to enter the alternative screen with the [`alternative`](@ref) function:

```julia-repl
julia> Screen.alternative() do
		   println("No scrollback!")
		   read(stdin, UInt8)
	   end
```

Ansillary will set raw mode (also known as [non-canonical mode](https://www.gnu.org/software/libc/manual/html_node/Canonical-or-Not.html)) when [`alternative`](@ref) is called as this is almost always what is wanted, to avoid this use the non-exported [`alternative!`](@ref) and [`standard!`](@ref) functions.

It also possible to only set raw mode using the [`raw`](@ref) function. The main benefits of raw mode are that the input stream is not line buffered allowing one byte to be read at a time (which in turn allows the application to respond to key presses), and that the input is not printed directly to the output allowing the application to handle printable input in it's own way (so that it can implement vim-style keybindings, for example).

To clear the screen Ansillary provides the [`clear!`](@ref) function, as well as the [`Area`](@ref) types for specifying which parts of the screen need clearing.

This short script

```julia
print("Some line...")
move!(Left(3))
clear!(FromCursorForward())
print("!")
```

will result in `Some line!` being printed to the terminal.

This module also provides the [`size`](@ref) as a slightly nicer wrapper around `displaysize`.
"""
module Screen

import ..Cursor
import ..TERMINAL

using REPL.Terminals: CSI, raw!

export All,
    CurrentLine,
    FromCursorBack,
    FromCursorDown,
    FromCursorForward,
    FromCursorUp,
    alternative,
    clear!,
    raw,
    size

"""
Permanently activate the alternative screen.

Use [`standard!`](@ref) to switch back to the standard screen.

This function **does not** set the terminal to raw mode, so that must be done manually where necessary.

!!! note

	You should prefer using [`alternative`](@ref) where possible as it's very easy to accidently leave the terminal on the alternative screen using this method.
"""
alternative!(terminal = TERMINAL[]) = print(terminal.out_stream, CSI, "?1049h")

"""
Permanently activate the standard screen.
"""
standard!(terminal = TERMINAL[]) = print(terminal.out_stream, CSI, "?1049l")

"""
Temporarily activate the alernative screen for the duration of the function.

This function also sets the terminal to raw mode as it is rare that you'll need the alternative screen but not raw mode. If the alternative screen is needed without setting raw mode, use [`alternative!`](@ref) and [`standard!`](@ref) directly.
"""
function alternative(f, terminal = TERMINAL[])
    raw(terminal) do
        alternative!(terminal)
        try
            f()
        finally
            standard!(terminal)
        end
    end
end

_clear!(terminal, code) = print(terminal.out_stream, CSI, code)

"""
An area of the screen.

See the documentation on it's subtypes for more details:

* [`All`](@ref)
* [`CurrentLine`](@ref)
* [`FromCursorBack`](@ref)
* [`FromCursorForward`](@ref)
* [`FromCursorDown`](@ref)
* [`FromCursorUp`](@ref)
"""
abstract type Area end

"""
Clear an area of the screen.

See subtypes of [`Area`](@ref) for more details.
"""
clear!(area::Area) = clear!(TERMINAL[], area)

"""
Clear the entire screen.
"""
struct All <: Area end
clear!(terminal, ::All) = _clear!(terminal, "2J")

"""
Clear the line that the cursor is currently on.
"""
struct CurrentLine <: Area end
clear!(terminal, ::CurrentLine) = _clear!(terminal, "2K")

"""
Clear from the start of the current line up to, and including, the cursor.
"""
struct FromCursorBack <: Area end
clear!(terminal, ::FromCursorBack) = _clear!(terminal, "1K")

"""
Clear from the cursor up to the end of the line.
"""
struct FromCursorForward <: Area end
clear!(terminal, ::FromCursorForward) = _clear!(terminal, "K")

"""
Clear from the cursor up to the end of the line and all lines below the cursor.
"""
struct FromCursorDown <: Area end
clear!(terminal, ::FromCursorDown) = _clear!(terminal, "J")

"""
Clear from the start of the current line up to the cursor and any lines above the cursor.
"""
struct FromCursorUp <: Area end
clear!(terminal, ::FromCursorUp) = _clear!(terminal, "1J")

"""
Temporarily set the terminal to raw mode for the duration of the function.

If switching to raw mode permanently is required use `REPL.Terminals.raw!`.
"""
function raw(f, terminal = TERMINAL[])
    raw!(terminal, true)
    try
        f()
    finally
        raw!(terminal, false)
    end
end

struct Size
    rows::UInt32
    columns::UInt32
end

"""
Get the current size of the terminal.

See also: `displaysize`.
"""
function size(terminal = TERMINAL[])
    (rows, columns) = displaysize(terminal.out_stream)
    Size(rows, columns)
end

end # module


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

using Markdown
using InteractiveUtils

# ╔═╡ 846cc37e-0a35-11eb-24ed-7f513230036a
using Pkg, DrWatson

# ╔═╡ 846cf9b6-0a35-11eb-3d6a-91ca6fd4e0d2
begin
	@quickactivate "StatisticalRethinkingTuring"
	using Turing
	using StatisticalRethinking
end

# ╔═╡ 8067a2c6-0a35-11eb-14cf-074eb476bf71
md"## Clip-05-45-49.jl"

# ╔═╡ 846d894e-0a35-11eb-365c-b95a0dc20b9f
md"### snippet 5.45"

# ╔═╡ 847e685e-0a35-11eb-3a5d-91ac3078947e
df = CSV.read(sr_datadir("Howell1.csv"), DataFrame);

# ╔═╡ 847ee568-0a35-11eb-1a4f-79f62939c240
md"### snippet 5.46"

# ╔═╡ 848e8644-0a35-11eb-2963-0d76c7c22c0b
begin
	μ_female = rand(Normal(178, 20), 10_000)
	μ_male = rand(Normal(178, 20), 10_000) .+ rand(Normal(0, 10), 10_000)
	Text(precis(DataFrame(; μ_female, μ_male); io=String))
end

# ╔═╡ 8496ffe0-0a35-11eb-2c3e-4159fcafa062
md"### snippet 5.47"

# ╔═╡ 849b6e0e-0a35-11eb-1451-279722a36141
df.sex = ifelse.(df.male .== 1, 2, 1);

# ╔═╡ 5da0339a-0a36-11eb-22d1-d9b7aea3644d
md"### snippet 5.48"

# ╔═╡ 6e0203fa-0a36-11eb-0eaa-357f2828ee4f
@model function m5_8(sex, height)
    σ ~ Uniform(0, 50)
    a ~ filldist(Normal(178, 20), 2)
    μ = a[sex]
    height ~ MvNormal(μ, σ)
end

# ╔═╡ 6e0238ac-0a36-11eb-2aad-316bea996d20
begin
	m5_8t = m5_8(df.sex, df.height)
	quap5_8t = quap(m5_8t, NelderMead())
	Text(precis(quap5_8t))
end

# ╔═╡ 6e02e748-0a36-11eb-313f-19358a19a2e1
md"### snippet 5.49"

# ╔═╡ 6e14cf6c-0a36-11eb-1b7c-61ee035c7e1d
begin
	post5_8t = DataFrame(rand(quap5_8t.distr, 1000)', quap5_8t.params)
	post5_8t.diff_fm = post5_8t[:, "a[1]"] .- post5_8t[:, "a[2]"]
	Text(precis(post5_8t; io=String))
end

# ╔═╡ c3a34dee-0b00-11eb-3225-6d90a2c31f7b
md"## End of clip-05-44-49t.jl"

# ╔═╡ Cell order:
# ╟─8067a2c6-0a35-11eb-14cf-074eb476bf71
# ╠═846cc37e-0a35-11eb-24ed-7f513230036a
# ╠═846cf9b6-0a35-11eb-3d6a-91ca6fd4e0d2
# ╟─846d894e-0a35-11eb-365c-b95a0dc20b9f
# ╠═847e685e-0a35-11eb-3a5d-91ac3078947e
# ╟─847ee568-0a35-11eb-1a4f-79f62939c240
# ╠═848e8644-0a35-11eb-2963-0d76c7c22c0b
# ╟─8496ffe0-0a35-11eb-2c3e-4159fcafa062
# ╠═849b6e0e-0a35-11eb-1451-279722a36141
# ╟─5da0339a-0a36-11eb-22d1-d9b7aea3644d
# ╠═6e0203fa-0a36-11eb-0eaa-357f2828ee4f
# ╠═6e0238ac-0a36-11eb-2aad-316bea996d20
# ╟─6e02e748-0a36-11eb-313f-19358a19a2e1
# ╠═6e14cf6c-0a36-11eb-1b7c-61ee035c7e1d
# ╟─c3a34dee-0b00-11eb-3225-6d90a2c31f7b


text:
module_path = splitdir(@__FILE__)[1]

# example 3 - static wind farm with n turbines

# NOTE: In Julia v0.7, JLD throws the error "Cannot `convert` an object of type
# getfield(JLD, Symbol("##JLD.AssociativeWrapper" unless, imported outside of
# WFVisual through the following line. I'm not sure, but I think it's because
# JLD must be imported before GeometricTools to be able to read Grid JLDs.
# Also, notice that `generate_windfarm(...)` takes an abnormal long time to run
# in v0.7.

using JLD
using LinearAlgebra
# Load WFVisual module
import WFVisual
using Statistics
wfv=WFVisual

# Load GeometricTools: https://github.com/byuflowlab/GeometricTools.jl
import GeometricTools
gt=GeometricTools

# Data path with geometry JLDs
#if Int(VERSION.major)==0 && Int(VERSION.minor)==6   # Case of Julia v0.6
 #   data_path = joinpath(module_path, "../data/")
#else                                                # Case of Julia v0.7
data_path = joinpath(module_path, "../datav07/")
#end

# Create save path where to store vtks
save_path = "example_3/"         # Save path of this example
gt.create_path(save_path, true)

# --------------------- WIND FARM LAYOUT ---------------------------------------
# turbine_x = [  0.0, 1000.0, -1000.0 ]

# turbine_y = [  0.0, 0.0, 0.0]
# turbine_z = [ 0.,  0.,  0.]

# hub_height = [ 100.,  100.,  100.]
# rotor_diameter = [ 100.,  100.,  100.]
# nBlades = [3, 3, 3]

# nturbs = length(turbine_x)
# hub_height = zeros(nturbs) .+ 100.0
# rotor_diameter = zeros(nturbs) .+ 100.0
# nBlades = zeros(Int8,nturbs) .+ 3

# wind_direction = 228.0

# yaw = [ 0.,  0.,  0.] .- wind_direction

turbine_x = [  909.98349606,  1005.0174903 ,   900.40238835,  -479.57607866,
                 937.92551703,   461.20472344,   123.06165965, -1073.3529325 ,
                -698.76200523, -1083.53094471,  1365.9338773 ,   204.10557045,
                 306.99725222,  -585.35767689,    10.9706543 ,   543.54927959,
                 -20.50212678, -1313.07762466, -1327.00852816,   137.06316521,
                1378.99410072,  -645.66250771,  -535.7876181 ,  -100.66538419,
               -1331.29278587]

turbine_y = [  500.59889721,  -974.65186327,    76.51586507,  1315.29789208,
                1039.37304649, -1321.85187113,  -300.57817461,  -765.82650713,
               -1213.14956771,   886.54973742,  -306.96866239,   194.51410928,
                 735.69027357,    81.95252777,  -892.14222889,  -930.41798152,
                 970.9117144 ,   398.75277565,    -3.52236011,  1393.27798468,
                 241.61586052,  -791.92492366,   919.26560593, -1396.37638078,
                -433.19849727]
turbine_z = [ 0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,
                0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.]

# hub_height = [ 100.,  100.,  100.,  100.,  100.,  100.,  100.,  100.,  100.,
#                 100.,  100.,  100.,  100.,  100.,  100.,  100.,  100.,  100.,
#                 100.,  100.,  100.,  100.,  100.,  100.,  100.]

# rotor_diameter = [ 100.,  100.,  100.,  100.,  100.,  100.,  100.,  100.,  100.,
#                 100.,  100.,  100.,  100.,  100.,  100.,  100.,  100.,  100.,
#                 100.,  100.,  100.,  100.,  100.,  100.,  100.]

# nBlades = [3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3,
#                3, 3]
nturbs = length(turbine_x)
hub_height = zeros(nturbs) .+ 100.0
rotor_diameter = zeros(nturbs) .+ 100.0
nBlades = zeros(Int8,nturbs) .+ 3

wind_direction = 228.0

yaw = [ 0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,
        0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.] .- wind_direction

# --------------------- PERIMETER AND FLUID DOMAIN -----------------------------
NDIVSx = 20              # Cells in the parametric x-direction
NDIVSy = 20            # Cells in the parametric y-direction
NDIVSz = 10              # Cells in the geometric z-direction

# Dummy perimeter
Rper = 1500.0
perimeter_points = Rper.*[ [cos(a), sin(a), 0] for a in range(0, 2*pi, 179)]

# --------------------- GENERATE FLOW FIELD ------------------------------------

_zmin = 0
_zmax = maximum(hub_height) + 1.25*maximum(rotor_diameter)/2
z_off = 0.0
spl_s=0.01
fdom = wfv.generate_perimetergrid(perimeter_points,
                                  NDIVSx, NDIVSy, NDIVSz;
                                  z_min=_zmin+z_off, z_max=_zmax+z_off,
                                  verify_spline=false,
                                  spl_s=spl_s,
                                  save_path=nothing
                                )
# Guassian wake function
magVinf = 12.0                          # (m/s) wind speed
Vinf_Oaxis = gt.rotation_matrix2(0, 0, wind_direction) # Free-stream coordinate system
invVinf_Oaxis = Vinf_Oaxis'             # Inverse transformation
refH = mean(hub_height)                 # (m) reference height for shear layer

function wake(X)

  k = 0.0325 # we can let this vary
  CT = 8.0/9.0
  loss = 0
  alpha = 0.12

  # Calculate loss in freestream direction
  for i = 1:length(turbine_x)     # Iterate over turbines

      # Distance from this turbine in global coordinates
      DX = X - [turbine_x[i], turbine_y[i], turbine_z[i]+hub_height[i]]

      # Distance in free-stream coordinates
      dX = invVinf_Oaxis*DX
      dx = dX[1]

      if dx > 0.
          dy = dX[2]
          dz = dX[3]
          sigma = k*dx + rotor_diameter[i]/sqrt(8.0)
          this_loss = (1.0-sqrt(1.0-CT/(8*sigma^2/(rotor_diameter[i]^2))))*
                          exp(-1.0/2.0*(dy/sigma)^2)*exp(-1.0/2.0*(dz/sigma)^2)
          loss += this_loss^2
      end
  end
  loss = sqrt(loss)

  windspeed= magVinf
  # Uncomment this line adds the shear layer
  windspeed *= (X[3]/refH)^alpha

  # Return velocity vector in global coordinate system
  return Vinf_Oaxis*[ windspeed*(1-loss), 0, 0 ]
end

gt.calculate_field(fdom, wake, "wake", "vector", "node")

# --------------------- GENERATE WIND FARM -------------------------------------
wfv.generate_windfarm(rotor_diameter, hub_height, nBlades,
                              turbine_x, turbine_y, turbine_z,
                              yaw,
                              perimeter_points, fdom;
                              NDIVSx=NDIVSx, NDIVSy=NDIVSy, NDIVSz=NDIVSz,
                              save_path=save_path, spl_s=0.01,
                              data_path=data_path, paraview=false);
nothing


text:
##

using IncrementalInference
using Test

##

@testset "test endless cycle case, issue #754" begin

##

# testgraph from issue #754
fg = generateGraph_LineStep(5;
                            poseEvery=1,
                            landmarkEvery=5,
                            posePriorsAt=[0,2],
                            sightDistance=4)
#
getSolverParams(fg).graphinit = false
getSolverParams(fg).treeinit = true
getSolverParams(fg).useMsgLikelihoods = true

##

# tree = buildTreeReset!(fg)
# drawTree(tree, show=true)

##

# ENV["JULIA_DEBUG"] = :csm_4

smtasks = Task[]
hist = IIF.solveTree!(fg; smtasks=smtasks); #, recordcliqs=ls(fg));

##

# hists = fetchCliqHistoryAll!(smtasks);
# printCSMHistorySequential(hists)
# printCSMHistoryLogical(hists)

# ##

# fnc_ = repeatCSMStep!(hists, 4, 6)

##

for var in sortDFG(ls(fg))
  sppe = getVariable(fg,var) |> getPPE |> IIF.getPPESuggested
  println("Testing ", var,": ", sppe)
  @test isapprox(sppe[1], parse(Int,string(var)[end]), atol=0.2)
end

##

# linear octo
N=8
fg = generateGraph_LineStep(N;
                            graphinit=false,
                            poseEvery=1,
                            landmarkEvery=N+1,
                            posePriorsAt=[0],
                            landmarkPriorsAt=[],
                            sightDistance=N+1)

deleteFactor!.(fg, [Symbol("x$(i)lm0f1") for i=1:(N-1)])

getSolverParams(fg).graphinit = false
getSolverParams(fg).treeinit = true
getSolverParams(fg).useMsgLikelihoods = true

smtasks = Task[]
tree = IIF.solveTree!(fg; smtasks=smtasks);

##

# hists = fetchCliqHistoryAll!(smtasks);
# printCSMHistorySequential(hists, 2=>1:50)
# printCSMHistoryLogical(hists)

# ##

# csmc_ = repeatCSMStep!(hists, 2, 14)

##

for var in sortDFG(ls(fg))
    sppe = getVariable(fg,var) |> getPPE |> IIF.getPPESuggested
    println("Testing ", var,": ", sppe)
    @test isapprox(sppe[1], parse(Int,string(var)[end]), atol=0.4)
end

##

# Larger graph
fg = generateGraph_LineStep(15;
                            poseEvery=1,
                            landmarkEvery=3,
                            posePriorsAt=[0,7,12],
                            landmarkPriorsAt=[0,3],
                            sightDistance=2)

getSolverParams(fg).graphinit = false
getSolverParams(fg).treeinit = true
getSolverParams(fg).useMsgLikelihoods = true

smtasks = Task[]
tree = IIF.solveTree!(fg; smtasks=smtasks);

for var in sortDFG(ls(fg))
    sppe = getVariable(fg,var) |> getPPE |> IIF.getPPESuggested
    println("Testing ", var,": ", sppe)
    s = findfirst(r"\d", string(var))[1]
    @test isapprox(sppe[1], parse(Int,string(var)[s:end]), atol=0.45)
end

##

end


text:
@testset "Dynamics" begin
    T = 3
    nx = 2
    nu = 1
    nw = 0
    w_dim = [nw for t = 1:T]

    function pendulum(z, u, w)
        mass = 1.0
        lc = 1.0
        gravity = 9.81
        damping = 0.1
        [z[2], (u[1] / ((mass * lc * lc)) - gravity * sin(z[1]) / lc - damping * z[2] / (mass * lc * lc))]
    end

    function euler_implicit(y, x, u, w)
        h = 0.1
        y - (x + h * pendulum(y, u, w))
    end

    dt = Dynamics(euler_implicit, nx, nx, nu, nw=nw);
    dyn = [dt for t = 1:T-1]
    model = DynamicsModel(dyn, w_dim=w_dim)

    x1 = ones(nx)
    u1 = ones(nu)
    w1 = zeros(nw)
    X = [x1 for t = 1:T]
    U = [u1 for t = 1:T]
    W = [w1 for t = 1:T]
    idx_dyn = DirectTrajectoryOptimization.constraint_indices(dyn)
    idx_jac = DirectTrajectoryOptimization.jacobian_indices(dyn)
    dim_x = model.dim.x
    dim_u = model.dim.u
    d = zeros(DirectTrajectoryOptimization.num_con(dyn))
    j = zeros(DirectTrajectoryOptimization.num_jac(dyn))

    dt.val(dt.val_cache, x1, x1, u1, w1)
    @test norm(dt.val_cache - euler_implicit(x1, x1, u1, w1)) < 1.0e-8
    dt.jac(dt.jac_cache, x1, x1, u1, w1)
    jac_dense = zeros(dt.ny, dt.nx + dt.nu + dt.ny)
    for (i, ji) in enumerate(dt.jac_cache)
        jac_dense[dt.sp_jac[1][i], dt.sp_jac[2][i]] = ji
    end
    jac_fd = ForwardDiff.jacobian(a -> euler_implicit(a[nx + nu .+ (1:nx)], a[1:nx], a[nx .+ (1:nu)], w1), [x1; u1; x1])
    @test norm(jac_dense - jac_fd) < 1.0e-8

    DirectTrajectoryOptimization.eval_con!(d, idx_dyn, dyn, X, U, W)
    @test norm(vcat(d...) - vcat([euler_implicit(X[t+1], X[t], U[t], W[t]) for t = 1:T-1]...)) < 1.0e-8
    # info = @benchmark DirectTrajectoryOptimization.eval_con!($d, $idx_dyn, $dyn, $X, $U, $W)

    DirectTrajectoryOptimization.eval_jac!(j, idx_jac, dyn, X, U, W)
    s = DirectTrajectoryOptimization.sparsity_jacobian(dyn, dim_x, dim_u)
    jac_dense = zeros(DirectTrajectoryOptimization.num_con(dyn), DirectTrajectoryOptimization.num_xuy(dyn))
    for (i, ji) in enumerate(j)
        jac_dense[s[i][1], s[i][2]] = ji
    end

    @test norm(jac_dense - [jac_fd zeros(dyn[2].nx, dyn[2].nu + dyn[2].ny); zeros(dyn[2].ny, dyn[1].nx + dyn[1].nu) jac_fd]) < 1.0e-8
    # info = @benchmark DirectTrajectoryOptimization.eval_jac!($j, $idx_jac, $dyn, $X, $U, $W)

    idx_x = model.idx.x
    idx_u = model.idx.u
    idx_xu = model.idx.xu
    idx_xuy = model.idx.xuy

    nz = sum([t < T ? dyn[t].nx : dyn[t-1].ny for t = 1:T]) + sum([dyn[t].nu for t = 1:T-1])
    z = rand(nz)
    x = [zero(z[idx_x[t]]) for t = 1:T]
    u = [[zero(z[idx_u[t]]) for t = 1:T-1]..., zeros(0)]

    DirectTrajectoryOptimization.trajectory!(x, u, z, idx_x, idx_u)
    z̄ = zero(z)
    for (t, idx) in enumerate(idx_x)
        z̄[idx] .= x[t]
    end
    for (t, idx) in enumerate(idx_u)
        z̄[idx] .= u[t]
    end

    @test norm(z - z̄) < 1.0e-8
    # info = @benchmark DirectTrajectoryOptimization.trajectory!($x, $u, $z, $idx_x, $idx_u)
end


text:
###################################################
## geo/trips.jl
## types representing trips
###################################################

"""
    `GeoTrip` : data for one geographical trip
"""
struct GeoTrip
    "pick-up lonlat coordinates"
    pLon::Float32
    pLat::Float32

    "drop-off lonlat coordinates"
    dLon::Float32
    dLat::Float32

    "Pick-up date and time"
    pTime::DateTime

    "Trip time (seconds)"
    time::Float32
end
function Base.show(io::IO, t::GeoTrip)
    @printf(io, "%s -- (%.5f,%.5f)=>(%.5f,%.5f): %7.1fm in %4.0fs", t.pTime, t.pLon, t.pLat, t.dLon, t.dLat, tripDistance(t), t.time)
end

"""
    `GeoData`: list of `GeoTrip`
"""
const GeoData = Vector{GeoTrip}

"""
    `stats`: print information
"""
function stats(t::GeoData)
    println("GeoData: trip information in GPS coordinates")
    println("--- $(length(t)) trips: $(sizeof(t)/1e6)MB")
end


text:
using LinearAlgebra, StaticArrays

"""
implements inverse kinematics for delta robot
"""
function inv_kin(x::Vector{Real}, p::DeltaRobotParams)
    @assert length(x) == 3
end

"""
computes the forward kinematics of the delta robot from the measured joint parameters phi, an initial guess
of the joint parameters x0, and the delta robot parameters p; the user may optionally pass in 2 additional
keyword arguments, namely abstol, which determines the absolute tolerance of the algorithm, and the maximum
number of iterations one would wish to take place

ex:
ee1, ee2, ee3, h1, h2, i, eps = for_kin(phi, [0.01,0.01,0.01,0.01,0.01,0.01], p)

- note on usage: if the initial guess vector is the zero vector, then one of the rows of the jacobian
  matrix becomes 0, which makes the inversion ill-posed
    - always initialize the matrix with a "good guess" otherwise with at least 1 non-zero element
"""
function for_kin(phi::Array{T, 1}, x0::Array{T, 1}, p::DeltaRobotParams; abstol=1e-15, max_iter=75) where T <: Real
    @assert length(phi) == 3
    @assert length(x0) == 6

    q11 = phi[1]
    q21 = phi[2]
    q31 = phi[3]
    ϵ = nothing
    l1 = p.l1
    l2 = p.l2
    r_base = p.rb
    r_platform = p.rp
    i = 0

    while (ϵ === nothing || abs(ϵ) > abstol) && i < max_iter
        q12 = x0[1];
        q13 = x0[2];
        q22 = x0[3];
        q23 = x0[4];
        q32 = x0[5];
        q33 = x0[6];

        J = @SMatrix [
            l2*cos(q12)*sin(q13) l2*cos(q13)*sin(q12) (l2*(cos(q22)*sin(q23) + 3^(1/2)*sin(q22)))/2 (l2*cos(q23)*sin(q22))/2 0 0;
            -l2*sin(q12) 0 -(l2*(sin(q22) - 3^(1/2)*cos(q22)*sin(q23)))/2 (3^(1/2)*l2*cos(q23)*sin(q22))/2 0 0;
            l2*cos(q12)*cos(q13) -l2*sin(q12)*sin(q13) -l2*cos(q22)*cos(q23) l2*sin(q22)*sin(q23) 0 0;
            l2*cos(q12)*sin(q13) l2*cos(q13)*sin(q12) 0 0 (l2*(cos(q32)*sin(q33) - 3^(1/2)*sin(q32)))/2 (l2*cos(q33)*sin(q32))/2;
            -l2*sin(q12) 0 0 0 -(l2*(sin(q32) + 3^(1/2)*cos(q32)*sin(q33)))/2 -(3^(1/2)*l2*cos(q33)*sin(q32))/2;
            l2*cos(q12)*cos(q13) -l2*sin(q12)*sin(q13) 0 0 -l2*cos(q32)*cos(q33) l2*sin(q32)*sin(q33)
        ];

        h = SVector{6}([
            l2*sin(q12)*sin(q13) - (3^(1/2)*(r_base - r_platform + l1*cos(q21) + l2*cos(q22)))/2 + (l2*sin(q22)*sin(q23))/2;
            (3*r_base)/2 - (3*r_platform)/2 + l1*cos(q11) + l2*cos(q12) + (l1*cos(q21))/2 + (l2*cos(q22))/2 + (3^(1/2)*l2*sin(q22)*sin(q23))/2;
            l1*sin(q11) - l1*sin(q21) + l2*cos(q13)*sin(q12) - l2*cos(q23)*sin(q22);
            (3^(1/2)*(r_base - r_platform + l1*cos(q31) + l2*cos(q32)))/2 + l2*sin(q12)*sin(q13) + (l2*sin(q32)*sin(q33))/2;
            (3*r_base)/2 - (3*r_platform)/2 + l1*cos(q11) + l2*cos(q12) + (l1*cos(q31))/2 + (l2*cos(q32))/2 - (3^(1/2)*l2*sin(q32)*sin(q33))/2;
            l1*sin(q11) - l1*sin(q31) + l2*cos(q13)*sin(q12) - l2*cos(q33)*sin(q32)
        ]);

        xk = x0 - J \ h;

        ϵ = sum(xk - x0);
        x0 = xk;
        i += 1;
    end

    x0 = convert(Vector{Real}, x0)
    x0 = [q11, x0[1], x0[2], q21, x0[3], x0[4], q31, x0[5], x0[6]];

    return x0, i, ϵ;
end

function joint_angles_to_ee(x0::Vector{Real})
    @assert length(x) == 9

    ee1 = [l2 * sin(x0[2]) * sin(x0[3]); r_base + l1 * cos(x0[1]) + l2 * cos(x0[2]); l1 * sin(x0[1]) + l2 * sin(x0[2]) * cos(x0[3])];
    ee2 = [l2 * sin(x0[5]) * sin(x0[6]); r_base + l1 * cos(x0[4]) + l2 * cos(x0[5]); l1 * sin(x0[5]) + l2 * sin(x0[5]) * cos(x0[6])];
    ee3 = [l2 * sin(x0[8]) * sin(x0[9]); r_base + l1 * cos(x0[7]) + l2 * cos(x0[8]); l1 * sin(x0[7]) + l2 * sin(x0[8]) * cos(x0[9])];

    rot_z(theta) = [
        cos(theta) -sin(theta) 0;
        sin(theta) cos(theta) 0;
        0 0 1
    ];
    Rz_2pi3 = rot_z(2 * pi / 3);
    Rz_neg_2pi3 = rot_z(-2 * pi / 3);
    h1 = ee1 - [0;r_platform;0] - Rz_neg_2pi3 * (ee2 - [0;r_platform;0]);
    h2 = ee1 - [0;r_platform;0] -     Rz_2pi3 * (ee3 - [0;r_platform;0]);

    @assert isapprox(sum(h1), 0.0; atol=1e-12)
    @assert isapprox(sum(h2), 0.0; atol=1e-12)

    return ee1 - [0;r_platform;0]
end

max_stars_repo_name: zborffs/Delta

max_issues_repo_name: zborffs/Delta

max_forks_repo_name: zborffs/Delta


text:
using DataFrames
import CSV
using EcologicalNetworks
using StatsPlots
using LinearAlgebra
using Statistics

# Functions
include("lib.jl")

# Read clover
include("read_clover.jl")

# Get array parameters
imputation_rank = parse(Int64, get(ENV, "SLURM_ARRAY_TASK_ID", "1"))

# Get them matrices
A = Array(clover.edges)

# Pools
filled = findall(A)
empty = findall(.!(A))

# Positions to sample
n_samples = Threads.nthreads()*20

dfs = [DataFrame(value = Bool[], position = [], model = Symbol[], initial = Float64[], updated = Float64[]) for i in 1:Threads.nthreads()]

QUER = [Float64.(A) for i in 1:Threads.nthreads()]
TMPLn1 = [n1(Q) for Q in QUER]
TMPLn2 = [n2(QUER[i]) for i in 1:Threads.nthreads()]
TMPLn3 = [n3(QUER[i]) for i in 1:Threads.nthreads()]

empties = rand(empty, n_samples)
filleds = rand(filled, n_samples)

Threads.@threads for i in 1:n_samples
    Qu = QUER[Threads.threadid()]
    try
        thisone = impute(Qu, TMPLn1[Threads.threadid()], empties[i]; rank=imputation_rank)
        push!(dfs[Threads.threadid()], (false, i, :connectance, thisone.initial, thisone.updated))
        thisone = impute(Qu, TMPLn2[Threads.threadid()], empties[i]; rank=imputation_rank)
        push!(dfs[Threads.threadid()], (false, i, :degree, thisone.initial, thisone.updated))
        thisone = impute(QUER[Threads.threadid()], TMPLn3[Threads.threadid()], empties[i]; rank=imputation_rank)
        push!(dfs[Threads.threadid()], (false, i, :hybrid, thisone.initial, thisone.updated))

        thisone = impute(Qu, TMPLn1[Threads.threadid()], filleds[i]; rank=imputation_rank)
        push!(dfs[Threads.threadid()], (true, i, :connectance, thisone.initial, thisone.updated))
        thisone = impute(Qu, TMPLn2[Threads.threadid()], filleds[i]; rank=imputation_rank)
        push!(dfs[Threads.threadid()], (true, i, :degree, thisone.initial, thisone.updated))
        thisone = impute(Qu, TMPLn3[Threads.threadid()], filleds[i]; rank=imputation_rank)
        push!(dfs[Threads.threadid()], (true, i, :hybrid, thisone.initial, thisone.updated))
    catch
        continue
    end
    CSV.write("svd-tuning-thread-$(Threads.threadid())-rank-$(imputation_rank).csv", dfs[Threads.threadid()])
end

max_stars_repo_name: viralemergence/trefle

max_issues_repo_name: viralemergence/trefle

max_forks_repo_name: viralemergence/trefle