text:
using Statistics, FFTW, LinearAlgebra#, JLD
# using BSON: @save, @load
# using DelimitedFiles, ArgParse,
include("./inNout.jl")
include("./corrFunctions.jl")
include("./test.jl")

mutable struct structureFactor
	ωT_list
	freqT_list
	CT_mixedWindows
	CT_arr

	ωL_list
	freqL_list
	CL_mixedWindows
	CL_arr

	ω_list
	freq_list
	C_mixedWindows
	C_arr

	CL_Im_mixedWindows
	CL_Im_arr

	CT_Im_mixedWindows
	CT_Im_arr
end

#=
#Loading Functions
=#

# function arrange_arrays(data, N)
#     #=
#     This arranges the data as vectors of x,y,z,vx,vy and vz
#     =#
#     x_indices = [2 + 6*i for i in 0:N-1]
#     y_indices = [3 + 6*i for i in 0:N-1]
#     z_indices = [4 + 6*i for i in 0:N-1]
#     vx_indices = [5 + 6*i for i in 0:N-1]
#     vy_indices = [6 + 6*i for i in 0:N-1]
#     vz_indices = [7 + 6*i for i in 0:N-1]
#
#     x, y, z = data[:,x_indices], data[:,y_indices], data[:,z_indices]
#     vx, vy, vz = data[:,vx_indices], data[:,vy_indices], data[:,vz_indices]
#     x, y, z, vx, vy, vz
# end

# function arrange_arrays(idx::Int, N)
#     pos, vel = load_data(i=idx, N=N)
#
#     x = hcat([pos[:,1+3*i] for i in 0:N-1]...)
#     y = hcat([pos[:,2+3*i] for i in 0:N-1]...)
#     z = hcat([pos[:,3+3*i] for i in 0:N-1]...)
#     vx = hcat([vel[:,1+3*i] for i in 0:N-1]...)
#     vy = hcat([vel[:,2+3*i] for i in 0:N-1]...)
#     vz = hcat([vel[:,3+3*i] for i in 0:N-1]...)
#     x, y, z, vx, vy, vz
# end

# function load_data(; i=1, path=PATH, N=num_part)
#     cd(path)
#     files = readdir(path * "/$i")
#     positions = readdlm(path * "/$i/positions.dat")
#     velocities = readdlm(path * "/$i/velocity.dat")
#     Array(reshape(positions, 3N, :)'), Array(reshape(velocities, 3N, :)')
# end

#=
Tests
=#
function test_data(;Δt = 0.01, steps = 5000, thermal_steps=1000,
			N=300, L = 1, dims=3, arrange = "Crystal", typeHessian="random")

	r, v, M = generateData(N; steps, dims, Δt, arrange, typeHessian)
	((r[1,thermal_steps+1:end,:]/L, r[2,thermal_steps+1:end,:]/L, r[3,thermal_steps+1:end,:]/L,
	v[1,thermal_steps+1:end,:], v[2,thermal_steps+1:end,:], v[3,thermal_steps+1:end,:]), M)
end

function partition_data(d, t_sample, dict)
    t_max = size(d[1],1)
    samp = Int(floor(t_max/t_sample))
    data = [(d[1][r,:],d[2][r,:],d[3][r,:],d[4][r,:],d[5][r,:],d[6][r,:]) for r in Iterators.partition(1:t_max,t_sample)][1:samp]
    @info "Data partitioned"
	append!(dict[:partitions][:t_samp], t_sample)
	append!(dict[:partitions][:samples], samp)
	append!(dict[:partitions][:t_max], samp * t_sample)
    data, samp * t_sample
end

function get_Corr_over_all_samples(d, dict; in_Fourier=false, flag="Transversal")

    data = d[1]
	N = dict[:num_part]
	t_max = dict[:partitions][:t_max][end]
	t_samp = dict[:partitions][:t_samp][end]
	samp = dict[:partitions][:samples][end]
	kmin, kmax = dict[:knum]
	L = dict[:L_char]
	m = dict[:mass]

	@info "Computing $flag part"
	if flag == "Transversal"
		@info "Sample 1"
	    CT, CTIm = getTransversalCorrelationFunction(data, N, t_samp; kmin = kmin, kmax = kmax, L = L, m=m)
	    CT_av, CTIm_av = reshape(CT,1,kmax-kmin+1,t_samp), reshape(CTIm,1,kmax-kmin+1,t_samp)
	    for i in 2:samp
	        @info "Sample $i"
	        data = d[i]
	        CT, CTIm = getTransversalCorrelationFunction(data, N, t_samp; kmin = kmin, kmax = kmax, L = L, m=m)

	        CT_av = vcat(CT_av, reshape(CT,1,kmax-kmin+1,t_samp))
	        CTIm_av = vcat(CTIm_av, reshape(CTIm,1,kmax-kmin+1,t_samp))
	    end
		return do_the_fou(CT_av, CTIm_av, t_samp; in_Fourier, L)
	elseif flag == "Longitudinal"
		@info "Sample 1"
		CL, CLIm = getLongitudinalCorrelationFunction(data, N, t_samp; kmin = kmin, kmax = kmax, L = L, m=m)
	    CL_av, CLIm_av = reshape(CL,1,kmax-kmin+1,t_samp), reshape(CLIm,1,kmax-kmin+1,t_samp)
	    for i in 2:samp
	        @info "Sample $i"
	        data = d[i]
			CL, CLIm = getLongitudinalCorrelationFunction(data, N, t_samp; kmin = kmin, kmax = kmax, L = L, m=m)

			CL_av = vcat(CL_av, reshape(CL,1,kmax-kmin+1,t_samp))
	        CLIm_av = vcat(CLIm_av, reshape(CLIm,1,kmax-kmin+1,t_samp))
	    end
		return do_the_fou(CL_av, CLIm_av, t_samp; in_Fourier, L)
	end
	# return CT_av, CTIm_av
end

function do_the_fou(CT_av, CTIm_av, tmax; kmin = 1, kmax = 10,
			L = 1, m=1, in_Fourier=true)
	if in_Fourier
        CT_av, CTIm_av = reshape(mean(CT_av,dims=1),kmax-kmin+1,tmax), reshape(mean(CTIm_av,dims=1),kmax-kmin+1,tmax)
        CT_FT = abs.(hcat(fft.([CT_av[i,:] for i in 1:kmax-kmin+1])...))
        CTIm_FT = abs.(hcat(fft.([CTIm_av[i,:] for i in 1:kmax-kmin+1])...))
        # return CT_FT, CTIm_FT
		return CT_FT, Array(transpose(CTIm_av))
    else
        return reshape(mean(CT_av,dims=1),kmax-kmin+1,tmax), reshape(mean(CTIm_av,dims=1),kmax-kmin+1,tmax)
    end
end

function get_ω_sample(data, dict, c::Bool)

    ωT_list = zeros(size(dict[:t_samp_list],1),dict[:knum][2]-dict[:knum][1] + 1)
	ωL_list = zeros(size(dict[:t_samp_list],1),dict[:knum][2]-dict[:knum][1] + 1)
	ωTim_list = zeros(size(dict[:t_samp_list],1),dict[:knum][2]-dict[:knum][1] + 1)
	ωLim_list = zeros(size(dict[:t_samp_list],1),dict[:knum][2]-dict[:knum][1] + 1)
	ω_list = zeros(size(dict[:t_samp_list],1),dict[:knum][2]-dict[:knum][1] + 1)
	CT_arr = zeros(maximum(dict[:t_samp_list]),dict[:knum][2]-dict[:knum][1] + 1,2,size(dict[:t_samp_list],1))
	CL_arr = zeros(maximum(dict[:t_samp_list]),dict[:knum][2]-dict[:knum][1] + 1,2,size(dict[:t_samp_list],1))
	CT_Im_arr = zeros(maximum(dict[:t_samp_list]),dict[:knum][2]-dict[:knum][1] + 1,2,size(dict[:t_samp_list],1))
	CL_Im_arr = zeros(maximum(dict[:t_samp_list]),dict[:knum][2]-dict[:knum][1] + 1,2,size(dict[:t_samp_list],1))
	C_arr = zeros(maximum(dict[:t_samp_list]),dict[:knum][2]-dict[:knum][1] + 1,2,size(dict[:t_samp_list],1))

	for (j, ts) in enumerate(dict[:t_samp_list])
        # ts = t_samp_list[j]
        d_part, t_m = partition_data(data, ts, dict)
        @info "$j : t_samp set to $ts timesteps. $(dict[:partitions][:samples][end]) samples for given t_samp."
        if dict[:partitions][:t_samp][end] * dict[:partitions][:samples][end] != dict[:partitions][:t_max][end]
            @warn "t_sample x n_sample does NOT equal t_m"
            break
        end

        CT, CT_Im = get_Corr_over_all_samples(d_part, dict; in_Fourier = c, flag = "Transversal")
		CL, CL_Im = get_Corr_over_all_samples(d_part, dict; in_Fourier = c, flag = "Longitudinal")
		C_total = CT .+ CL

		ωlim = Int(floor(ts/2))

		gatherDataToArrays(j, ωlim, CT, ts, dict, ωT_list, CT_arr)
		gatherDataToArrays(j, ωlim, CL, ts, dict, ωL_list, CL_arr)
		gatherDataToArrays(j, ωlim, CT_Im, ts, dict, ωTim_list, CT_Im_arr, in_fourier=false)
		gatherDataToArrays(j, ωlim, CL_Im, ts, dict, ωLim_list, CL_Im_arr, in_fourier=false)
		gatherDataToArrays(j, ωlim, C_total, ts, dict, ω_list, C_arr)

    end
    freqT_list = hcat([(ωT_list[i,:] .- 1) ./(ts * dict[:dt]) for (i,ts) in enumerate(dict[:t_samp_list])]...)
	freqL_list = hcat([(ωL_list[i,:] .- 1) ./(ts * dict[:dt]) for (i,ts) in enumerate(dict[:t_samp_list])]...)
	freq_list = hcat([(ω_list[i,:] .- 1) ./(ts * dict[:dt]) for (i,ts) in enumerate(dict[:t_samp_list])]...)

	CT_mixedWindows = getMixedWindows(CT_arr, dict)
	CL_mixedWindows = getMixedWindows(CL_arr, dict)
	# CT_Im_mixedWindows = getMixedWindows(CT_Im_arr, dict)
	# CL_Im_mixedWindows = getMixedWindows(CL_Im_arr, dict)
	C_mixedWindows = getMixedWindows(C_arr, dict)

	# ωT_list, freqT_list, CT_mixedWindows, CT_arr, ωL_list, freqL_list, CL_mixedWindows, CL_arr
	structureFactor(ωT_list, freqT_list, CT_mixedWindows, CT_arr, ωL_list, freqL_list, CL_mixedWindows, CL_arr, ω_list,	freq_list, C_mixedWindows, C_arr, undef, CL_Im_arr, undef, CT_Im_arr)
end

function gatherDataToArrays(j, ωlim, C, ts, dict, ω_list, C_arr; in_fourier=true)
	if in_fourier
		ω_list[j,:] = vcat([findall(x->x == maximum(C[1:ωlim,i]), C[1:ωlim,i]) for i in 1:size(C,2)]...)
		freq = hcat([[l for l in 0:size(C,1)-1] ./ (ts * dict[:dt]) for i in 1:size(C,2)]...)
		C_arr[1:size(C,1),:,:,j] = cat(C,freq, dims=3)
	else
		C_arr[1:size(C,1),:,:,j] = cat(C, hcat([dict[:dt] .* collect(1:size(C,1)) for i in 1:size(C,2)]...), dims=3)
	end
end

function getMixedWindows(C_arr, dict)
	C_mixedWindows = arrange_Corr(C_arr; len=Int(floor(minimum(dict[:t_samp_list])/2)), Δk = dict[:knum][2]-dict[:knum][1] + 1)
	C_mixedWindows = envelope.([C_mixedWindows[i,:,:] for i in 1:dict[:knum][2]-dict[:knum][1] + 1])
	C_mixedWindows = permutedims(cat(C_mixedWindows...,dims=3),(3,1,2))
	C_mixedWindows
end

function arrange_Corr(CT_arr; len=30, Δk=10)
        t_len=size(CT_arr,4)
        l = [i*len for i in 1:t_len]
        new_CT = zeros(Δk,sum(l),2)
        for k in 1:Δk
                x = vcat([CT_arr[1:i*len,k,2,i] for i in 1:t_len]...)
                y = vcat([CT_arr[1:i*len,k,1,i] for i in 1:t_len]...)
                idx = sortperm(x)
                new_CT[k,:,:] = hcat([x[idx],y[idx]]...)
        end
        new_CT
end

function envelope(ar)
    arr = ar[:,1]
    arr2 = ar[:,2]

    l = zeros(size(arr))
    ind = 1
    for i in 1:size(l,1)
        if ind <= size(arr2,1)
            ind = findall(x->x == maximum(arr2[arr[ind] .== arr]), arr2)[1]
            l[i] = ind
            ind = maximum(indexin(arr2[arr[ind] .== arr], arr2)) + 1
        end
    end
    ll = Array{Int,1}(l[l .>0])
    ar[ll,:]
end

function get_freq_max(CT_arr)
        freq_max = zeros(10,2)
        for k in 1:10
                idx = sortperm(CT_arr[k,:,2])[end]
                freq_max[k,:] = [k, CT_arr[k,idx,1]]
        end
        freq_max
end

function get_freqs(dict; test=false, ω = 0.0002, Δt = 10.0, δ = 0.5)
    if test
		t_MAX = dict[:t_max]
	    L_char = dict[:L_char]
	    num_part = dict[:num_part]
		t_samp_list = dict[:t_samp_list]
		d, M = test_data(steps = Int(t_MAX/dict[:dt]), thermal_steps=0, L = L_char,
		        Δt = dict[:dt], N=num_part, arrange = "Crystal", typeHessian="random")
    else
        d = load_data(dict)
    end

    get_ω_sample(d, dict, true)
end

function main(dict; test=false)
	# dim = dict[:dim]
	@info "Starting..."
	# ω_list, freq_list, CT_mixedWindows, CT_arr, ωL_list, freqL_list, CL_mixedWindows, CL_arr = get_freqs(dict)
	s = get_freqs(dict; test=test)
	@info "Computation done"

	@info "Generating Plots..."
	savePlots(s.ωT_list, s.freqT_list, s.CT_mixedWindows, s.CT_arr, dict, flag="T")
	savePlots(s.ωL_list, s.freqL_list, s.CL_mixedWindows, s.CL_arr, dict, flag="L")
	savePlots(s.ω_list, s.freq_list, s.C_mixedWindows, s.C_arr, dict, flag="total")
	savePlots(s.CL_Im_arr, dict, flag="L_im")
	savePlots(s.CT_Im_arr, dict, flag="T_im")

	@info "Saving Data..."
	saveData(s.ωT_list, s.freqT_list, s.CT_mixedWindows, s.CT_arr, dict, flag="T")
	saveData(s.ωL_list, s.freqL_list, s.CL_mixedWindows, s.CL_arr, dict, flag="L")
	saveData(s.ω_list, s.freq_list, s.C_mixedWindows, s.C_arr, dict, flag="total")
	saveData(s.CL_Im_arr, dict, flag="L_im")
	saveData(s.CT_Im_arr, dict, flag="T_im")
	@info "Done!"
	writedlm(dict[:PathOut] * "/" * split(dict[:title],".")[1] * "_freq_mixedWindows_T.dat",
			get_freq_max(s.CT_mixedWindows))
	writedlm(dict[:PathOut] * "/" * split(dict[:title],".")[1] * "_freq_mixedWindows_L.dat",
			get_freq_max(s.CL_mixedWindows))
	writedlm(dict[:PathOut] * "/" * split(dict[:title],".")[1] * "_freq_mixedWindows_total.dat",
			get_freq_max(s.C_mixedWindows))
end

max_stars_repo_name: jquetzalcoatl/DynamicStructureFactor

max_issues_repo_name: jquetzalcoatl/DynamicStructureFactor

max_forks_repo_name: jquetzalcoatl/DynamicStructureFactor


text:
# Create example data
x = LinRange(0, 3, 20)
y = [sin.(x) exp.(-x)]
# Set a new color scheme and plot the data
subplot(1,2,1)
colorscheme("light")
plot(x, y)
# Make a second plot with a particular scheme
subplot(1,2,2)
plot(x, y, scheme=3) # solarized light
# Now change the global scheme and redraw
# (this only affects the first plot)
colorscheme("solarized dark")
draw(gcf())

ext: jl


text:
function T=TFromThetaRho(Th,Rho,Param)
T=(Param.Rd*Th/Param.p0^Param.kappa).^(1/(1-Param.kappa))./(Param.Rd*Rho);
end


text:
julia_version = v"1.3.1"
include("../common.jl")

max_stars_repo_name: waralex/Yggdrasil

max_issues_repo_name: waralex/Yggdrasil

max_forks_repo_name: waralex/Yggdrasil


text:
struct ContactStructure
	distribution # store the hyperplane dist. as R^3 -> R^6 function
end

# the standard contact structure (x_1 dx_2 + dx_3)
standard_structure = ContactStructure((x,y,z) -> [0 1 -x;1 0 0])

# @TODO: Turn into something nicer
function structuralize_form(α::OneForm)
	ContactStructure((x,y,z) -> extract_from_oneform(α,x,y,z))
end

function extract_from_oneform(α::OneForm,x,y,z, tol = 1e-4)
	mat = [α.rep(x,y,z); α.rep(x,y,z); α.rep(x,y,z)]
	basis = kernel_basis(mat)
	if length(basis) != 2
		throw_error = true
		# check if the vectors are not actually the same solution twice
		if length(basis) == 3
			for i ∈ [1,2,3]
				for j ∈ setdiff([1,2,3],i)
					if sum(abs.(basis[i]-basis[j])) < tol || sum(abs.(basis[i]+basis[j])) < tol
						throw_error = false
						# Make sure that the two unique vectors are the output
						basis[1], basis[2] = basis[i], basis[setdiff([1,2,3],[i,j])[1]]
						break
					end
				end
			end
		end
		if throw_error
			throw(ArgumentError("'Contact form' has kernel of dimension $(length(basis)) which is wrong. Must be 2."))
		else
			return [transpose(basis[1]);transpose(basis[2])]
		end
	else
		return [transpose(basis[1]);transpose(basis[2])]
	end
end


text:
abstract type AbstractFoo end
abstract type AbstractPerson end
abstract type AbstractStudent <: AbstractPerson end
abstract type AbstractStudent2 <: AbstractPerson end
mutable struct Foo <: AbstractFoo
end
function bar(self::Foo)::Int64
    return baz(self)
end

function baz(self::Foo)::Int64
    return 10
end

function bar_str(self::Foo)::String
    return "a"
end

mutable struct Person <: AbstractPerson
    name::String
end
function get_name(self::Person)::String
    return self.name
end

mutable struct Student <: AbstractStudent
    name::String
    student_number::Int64
    domain::String

    Student(name::String, student_number::Int64, domain::String = "school.student.pt") =
        new(name, student_number, domain)
end
function get_name(self::Student)
    return "$(self.name) - $(self.student_number)"
end

mutable struct Student2 <: AbstractStudent2
    name::String
    student_number::Int64
    domain::String

    Student2(name::String, student_number::Int64, domain::String = "school.student.pt") =
        begin
            if student_number < 0
                throw(ValueError("Student number must be a positive number"))
            end
            new(name, student_number, domain)
        end
end

if abspath(PROGRAM_FILE) == @__FILE__
    f = Foo()
    b = bar(f)
    @assert(b == 10)
    c = bar_str(f)
    @assert(c == "a")
    p = Person("P")
    s = Student("S", 111111)
    @assert(get_name(p) == "P")
    @assert(get_name(s) == "S - 111111")
    s2 = Student2("S2", 123)
    println("OK")
end


text:
# Template providing lexer runtime functions

# Parses str with table from index i
# Returns (accepted, end_pos, next_pos)
function _lparse(table, str, pos)
    state = 1
    last_pos = deepcopy(pos)
    while pos.index <= length(str)
        char = str[pos.index]
        tok = table.tok2index(char)

        nstate = table.actions[state, tok]
        if nstate == nothing
            break
        end

        state = nstate
        last_pos = deepcopy(pos)

        # Change pos
        pos.index += 1
        if char == '\n'
            pos.line += 1
            pos.column = 1
        else
            pos.column += 1
        end
    end

    return (table.terminal_states[state], last_pos, pos)
end

# Lexer parse
# regs : Tuple(id, table, rule)
# Returns all tokens
function lparse(src, regs)
    tokens = []
    pos = pos_new()
    while pos.index <= length(src)
        # Test each reg
        err = true
        newpos = nothing
        for (id, reg, rule) in regs
            # endpos is just before newpos
            acc, endpos, newpos = _lparse(reg, src, deepcopy(pos))

            # Match
            if acc && newpos != pos
                err = false

                # Send token if not ignored
                if rule != nothing
                    data = rule(src[pos.index:newpos.index - 1])
                    tok = tok_new(id, terminal=true, data=data, start_pos=pos, end_pos=endpos)

                    push!(tokens, tok)
                end
                break
            end
        end

        # Error case (no match)
        if err
            line_start = findprev('\n', src, pos.index)
            line_end = findnext('\n', src, pos.index + 1)

            line_start == nothing && (line_start = 0)
            line_end == nothing && (line_end = length(src) + 1)

            println("--- Error Info ---")
            if line_start < line_end
                println(src[line_start + 1:line_end - 1])
                println(repeat(' ', pos.column - 1) * "^")
            end

            error("Invalid syntax at position $pos, no token matched")
        end

        pos = newpos
    end

    # End token
    push!(tokens, tok_end(start_pos=pos, end_pos=pos))

    return tokens
end


text:
__precompile__()
module GadgetGalaxies

using LinearAlgebra: rotate!

using Cosmology
using GadgetIO
using LinearAlgebra
using Statistics
using Unitful
using UnitfulAstro

include("snapshot.jl")
include("galaxy.jl")
include("read_halo.jl")
include("transformations.jl")
include("shapes.jl")
include("global_parameters.jl")
include("kinematic_parameters.jl")
include("units.jl")
include("utils.jl")

# structs
export Snapshot,
       Galaxy,
       GalaxyGroup,

       particle_type_id,

       # read halo
       read_halo!,
       read_redshift,
       read_header_particle_mass,
       read_galaxy_prop,
       read_galaxy_pos,
       read_galaxy_vel,
       is_main_halo,
       get_group,
       get_first_subhalo,

       # transformations
       rotate,
       rotate!,
       rotate_edgeon,
       rotate_edgeon!,
       translate,
       translate!,
       center_of_mass_iterative,
       translate_to_center_of_mass_iterative,
       translate_to_center_of_mass_iterative!,
       center_of_velocity,
       translate_to_center_of_velocity,
       translate_to_center_of_velocity!,

       # shapes
       Ellipsoid,
       Sphere,
       rotation_matrix_edgeon,
       rotation_matrix_2D,
       rotation_matrix_axis_ratios,
       r²_sphere,
       r²_circle,
       r²_ellipsoid,
       triaxiality,
       ellipticity,
       eccentricity,

       # global parameters
       half_mass_radius,
       half_mass_radius_2D,

       # kinematic parameters
       angular_momentum,
       b_value,
       circular_velocity,
       specific_angular_momentum,
       spin_parameter,

       # units
       convert_units!,
       convert_units_full,
       convert_units_physical,
       convert_units_physical!,

       simulation_units_pos,
       simulation_units_pos!,
       convert_units_pos,
       convert_units_full_pos,
       convert_units_physical_pos,
       convert_units_physical_pos!,
       simulation_units_vel,
       simulation_units_vel!,
       convert_units_vel,
       convert_units_full_vel,
       convert_units_physical_vel,
       convert_units_physical_vel!,
       simulation_units_temp,
       simulation_units_temp!,
       convert_units_temp,
       convert_units_full_temp,
       convert_units_physical_temp,
       convert_units_physical_temp!,
       simulation_units_mass,
       simulation_units_mass!,
       convert_units_mass,
       convert_units_full_mass,
       convert_units_physical_mass,
       convert_units_physical_mass!,
       simulation_units_age,
       simulation_units_age!,
       convert_units_age,
       convert_units_full_age,
       convert_units_physical_age,
       convert_units_physical_age!,

       convert_units_solar_metallicity

end


text:
# see tutorial https://turinglang.github.io/TuringCallbacks.jl/dev/

module TuringCLIExample

using TuringCallbacks: join
using Turing, TuringCallbacks, TensorBoardLogger, ArgParse, StatsPlots, CSV, DataFrames
import ArgParse.parse_item

function ArgParse.parse_item(::Type{Vector{Float64}}, x::AbstractString)
    return parse.(Float64, split(strip(x, [' ','\"',']','[']), ','))
end

function simulate_and_estimate(d)
    s_prior_alpha, s_prior_theta = d.prior
    @model function demo(x; s_prior_alpha, s_prior_theta)
        s ~ InverseGamma(s_prior_alpha, s_prior_theta)
        m ~ Normal(0, √s)
        for i in eachindex(x)
            x[i] ~ Normal(m, √s)
        end
    end

    xs = randn(100) .+ 1
    model = demo(xs; s_prior_alpha, s_prior_theta)

    # Sampling
    println("Generating $(d.num_samples) samples")
    callback = TensorBoardCallback(joinpath(pkgdir(TuringCLIExample),"tensorboard_logs/run"))
    alg = NUTS(d.num_adapts, d.target_acceptance_rate)
    chain = sample(model, alg, d.num_samples; callback)

    println("Generating trace plot")
    trace_plot = plot(chain, seriestype=:traceplot)
    savefig(trace_plot,  joinpath(callback.logger.logdir, "traceplots.png"))

    println("Summarizing the chain")
    sum_stats = describe(chain)
    param_names = sum_stats[1][:,1]
    param_mean = sum_stats[1][:,2]
    param_sd = sum_stats[1][:,3]
    param_ess = sum_stats[1][:,6]
    param_rhat = sum_stats[1][:,7]
    param_ess_per_sec = sum_stats[1][:, 8]

    CSV.write(
        joinpath(callback.logger.logdir, "summary.csv"),
        DataFrame(
            parameter=param_names,
            mean=param_mean,
            sd=param_sd,
            ess=param_ess,
            rhat=param_rhat,
            ess_per_sec= param_ess_per_sec
        )
    )

    # Log the ESS/sec and rhat.  Nice to show as summary results from tensorboard
    for (i, name) = enumerate(param_names)
        TensorBoardLogger.log_value(
            callback.logger,
            "$(name)_ess_per_sec",
            param_ess_per_sec[i],
        )
        TensorBoardLogger.log_value(
            callback.logger,
            "$(name)_rhat",
            param_rhat[i],
        )
    end
end

# Entry for script
function main(args = ARGS)
    d = parse_commandline(args)
    simulate_and_estimate((;d...)) # to named tuple
end

function parse_commandline(args)

    s = ArgParseSettings(fromfile_prefix_chars=['@'])

    # See src/defaults.txt
    @add_arg_table! s begin
        "--num_samples"
        help = "samples to draw in chain"
        arg_type = Int64
        "--num_adapts"
        help = "number of adaptations for NUTS"
        arg_type = Int64
        "--target_acceptance_rate"
        help = "Target acceptance rate for dual averaging."
        arg_type = Float64
        "--prior"
        help = "prior parameters in InverseGamma prior for s"
        arg_type = Vector{Float64}
    end

    args_with_default = vcat("@$(pkgdir(TuringCLIExample))/src/defaults.txt", args)
    return parse_args(args_with_default, s;as_symbols=true)
end
end #module

max_stars_repo_name: HighDimensionalEconLab/TuringCLIExample

max_issues_repo_name: HighDimensionalEconLab/TuringCLIExample

max_forks_repo_name: HighDimensionalEconLab/TuringCLIExample


text:
#push!(Base.LOAD_PATH,joinpath(dirname(@__FILE__),"..",".."))

using Base.Test, MathOptInterface
using LinQuadOptInterface

const MOIT = MathOptInterface.Test
const LQOI = LinQuadOptInterface

@testset "LinQuadOptInterface" begin
    @testset "Linear tests" begin
        linconfig = MOIT.TestConfig(solve = false)
        solver = LQOI.MockLinQuadOptimizer()
        MOIT.contlineartest(solver , linconfig, ["linear10"])
    end

    # @testset "Quadratic tests" begin
    #     quadconfig = MOIT.TestConfig(solve=false)
    #     solver = LQOI.MockLinQuadOptimizer()
    #     MOIT.contquadratictest(solver, quadconfig)
    # end

    @testset "Linear Conic tests" begin
        linconfig = MOIT.TestConfig(solve=false)
        solver = LQOI.MockLinQuadOptimizer()
        MOIT.lintest(solver, linconfig)
    end

    @testset "Integer Linear tests" begin
        intconfig = MOIT.TestConfig(solve=false)
        solver = LQOI.MockLinQuadOptimizer()
        MOIT.intlineartest(solver, intconfig, ["int2"])
    end

    @testset "ModelLike tests" begin
        intconfig = MOIT.TestConfig(solve=false)
        solver = LQOI.MockLinQuadOptimizer()
        MOIT.validtest(solver)
        MOIT.emptytest(solver)
        solver2 = LQOI.MockLinQuadOptimizer()
        MOIT.copytest(solver,solver2)
    end
end
;


text:
module adm
export HPDinterval

function HPDinterval(x::Int)
    x + 5
end

end


text:
using LinQuadOptInterface
const LQOI = LinQuadOptInterface
const MOI  = LQOI.MOI

const SUPPORTED_OBJECTIVES = [
    LQOI.Linear
    LQOI.SinVar
]

const SUPPORTED_CONSTRAINTS = [
    (LQOI.Linear, LQOI.EQ),
    (LQOI.Linear, LQOI.LE),
    (LQOI.Linear, LQOI.GE),
    (LQOI.SinVar, LQOI.EQ),
    (LQOI.SinVar, LQOI.LE),
    (LQOI.SinVar, LQOI.GE),
    (LQOI.SinVar, LQOI.IV),
    (LQOI.SinVar, MOI.ZeroOne),
    (LQOI.SinVar, MOI.Integer),
    (LQOI.VecVar, MOI.Nonnegatives),
    (LQOI.VecVar, MOI.Nonpositives),
    (LQOI.VecVar, MOI.Zeros),
    (LQOI.VecLin, MOI.Nonnegatives),
    (LQOI.VecLin, MOI.Nonpositives),
    (LQOI.VecLin, MOI.Zeros)
]

mutable struct Optimizer <: LQOI.LinQuadOptimizer
    LQOI.@LinQuadOptimizerBase
    env::Env # Cplex environment
    Optimizer(::Void) = new()
end

function LQOI.LinearQuadraticModel(::Type{Optimizer},env)
    env = Env()
    return Model(env::Env)
end

function Optimizer(;kwargs...)
    env = Env()
    for (name,value) in kwargs
        set_param!(env, string(name), value)
    end
    model = Optimizer(nothing)
    model.env = env
    MOI.empty!(model)
    return model
end

LQOI.supported_constraints(::Optimizer) = SUPPORTED_CONSTRAINTS
LQOI.supported_objectives(::Optimizer)  = SUPPORTED_OBJECTIVES

LQOI.backend_type(model::Optimizer, ::MOI.EqualTo{Float64})     = Cchar('E')
LQOI.backend_type(model::Optimizer, ::MOI.LessThan{Float64})    = Cchar('L')
LQOI.backend_type(model::Optimizer, ::MOI.GreaterThan{Float64}) = Cchar('G')

LQOI.backend_type(model::Optimizer, ::MOI.Zeros)                = Cchar('B')
LQOI.backend_type(model::Optimizer, ::MOI.Nonpositives)         = Cchar('L')
LQOI.backend_type(model::Optimizer, ::MOI.Nonnegatives)         = Cchar('E')

# TODO - improve single type
function LQOI.change_variable_bounds!(model::Optimizer,
        columns::Vector{Int}, values::Vector{Float64}, senses::Vector{Cchar})

    c_api_chgbds(model.inner, ivec(columns), senses, values)
end

function LQOI.get_variable_lowerbound(model::Optimizer, column::Int)
    c_api_getlb(model.inner, Cint(column), Cint(column))[1]
end

function LQOI.get_variable_upperbound(model::Optimizer, column::Int)
    c_api_getub(model.inner, Cint(column), Cint(column))[1]
end

function LQOI.get_number_linear_constraints(model::Optimizer)
    c_api_getnumrows(model.inner)
end

function LQOI.add_linear_constraints!(model::Optimizer,
        A::LQOI.CSRMatrix{Float64}, sense::Vector{Cchar}, rhs::Vector{Float64})

    c_api_addrows(model.inner, ivec(A.row_pointers), ivec(A.columns),
                  A.coefficients, sense, rhs)
end

function LQOI.get_rhs(model::Optimizer, row::Int)
    rhs = Vector{Cdouble}(1)
    c_api_getrhs(model.inner, rhs, Cint(row), Cint(row))
    return rhs[1]
end

function LQOI.get_linear_constraint(model::Optimizer, row::Int)
    (nzcnt, rmatbeg, rmatind, rmatval) =
            c_api_getrows(model.inner, Cint(row), Cint(row))
    return rmatind[1:nzcnt], rmatval[1:nzcnt]
end

function LQOI.change_matrix_coefficient!(model::Optimizer,
                                         row::Int, col::Int, coef::Float64)

    c_api_chgcoef(model.inner, Cint(row), Cint(col), coef)
end

function LQOI.change_objective_coefficient!(model::Optimizer, col::Int,
                                            coef::Float64)

    c_api_chgobj(model.inner, [Cint(col)], [coef])
end

function LQOI.change_rhs_coefficient!(model::Optimizer, row::Int,
                                      coef::Float64)

    c_api_chgrhs(model.inner, [Cint(row)], [coef])
end

function LQOI.delete_linear_constraints!(model::Optimizer,
                                         first_row::Int, last_row::Int)

    c_api_delrows(model.inner, Cint(first_row), Cint(last_row))
end

function LQOI.change_variable_types!(model::Optimizer,
        columns::Vector{Int}, vtypes::Vector{Cchar})

    c_api_chgctype(model.inner, ivec(columns), vtypes)
end

function LQOI.change_linear_constraint_sense!(model::Optimizer,
        rows::Vector{Int}, senses::Vector{Cchar})

    c_api_chgsense(model.inner, ivec(rows), senses)
end

function LQOI.set_linear_objective!(model::Optimizer,
        columns::Vector{Int}, coefficients::Vector{Float64})

    n = num_var(model.inner)
    all_coefs = zeros(Float64, n)
    for (col, coef) in zip(columns, coefficients)
        all_coefs[col] += coef
    end
    c_api_chgobj(model.inner, Cint[1:n;], all_coefs)
end

function LQOI.change_objective_sense!(model::Optimizer, symbol)
    if symbol == :min
        c_api_chgobjsen(model.inner, Cint(1))
    else
        c_api_chgobjsen(model.inner, Cint(-1))
    end
end

function LQOI.get_linear_objective!(model::Optimizer, x)
    c_api_getobj(model.inner, x, Cint(1), c_api_getnumcols(model.inner))
end

function LQOI.get_objectivesense(model::Optimizer)
    s = c_api_getobjsen(model.inner)
    if s == 1
        return MOI.MinSense
    else
        return MOI.MaxSense
    end
end

function LQOI.get_number_variables(model::Optimizer)
    return c_api_getnumcols(model.inner)
end

function LQOI.add_variables!(model::Optimizer, N::Int)
    add_vars!(model.inner, zeros(Float64, N), fill(-Inf, N), fill(Inf, N))
end

function LQOI.delete_variables!(model::Optimizer,
                                first_col::Int, last_col::Int)

    c_api_delcols(model.inner, Cint(first_col), Cint(last_col))
end

function LQOI.solve_mip_problem!(model::Optimizer)
    model.inner.has_int = true
    LQOI.solve_linear_problem!(model)
end

function LQOI.solve_linear_problem!(model::Optimizer)
    optimize!(model.inner)
end

function LQOI.get_termination_status(model::Optimizer)
    stat = c_api_getstat(model.inner)

    if stat == 1 # CPX_STAT_OPTIMAL
        return MOI.Success
    elseif stat == 3 # CPX_STAT_INFEASIBLE
        return MOI.InfeasibleNoResult
    elseif stat == 4 # CPX_STAT_INForUNBD
        return MOI.InfeasibleOrUnbounded
    elseif stat == 2 # CPX_STAT_UNBOUNDED
        return MOI.UnboundedNoResult
    elseif stat in [12, 21, 22, 36] # CPX_STAT_*ABORT*_OBJ_LIM
        return MOI.ObjectiveLimit
    elseif stat in [10,34] # CPX_STAT_*ABORT_IT_LIM
        return MOI.IterationLimit
    elseif stat == 53 # CPX_STAT_CONFLICT_ABORT_NODE_LIM
        return MOI.NodeLimit
    elseif stat in [11, 25, 33, 39] # CPX_STAT_*ABORT*TIME_LIM
        return MOI.TimeLimit
    elseif stat == 5 # CPX_STAT_OPTIMAL_INFEAS
        return MOI.NumericalError

    # MIP STATUS
    elseif stat in [101, 102] # CPXMIP_OPTIMAL, CPXMIP_OPTIMAL_TOL
        return MOI.Success
    elseif stat == 103 # CPXMIP_INFEASIBLE
        return MOI.InfeasibleNoResult
    elseif stat == 119 # CPXMIP_INForUNBD
        return MOI.InfeasibleOrUnbounded
    elseif stat == 118 # CPXMIP_UNBOUNDED
        return MOI.UnboundedNoResult
    elseif stat == [105, 106] # CPXMIP_NODE_LIM*
        return MOI.NodeLimit
    elseif stat in [107, 108, 131, 132] # CPXMIP_*TIME_LIM*
        return MOI.TimeLimit
    else
        return MOI.OtherError
    end
end

function LQOI.get_primal_status(model::Optimizer)
    stat = c_api_solninfo(model.inner)[3]
    if stat == 1
        return MOI.FeasiblePoint
    else
        return MOI.UnknownResultStatus
    end
end

function LQOI.get_dual_status(model::Optimizer)
    if model.inner.has_int
        return MOI.UnknownResultStatus
    end
    stat = c_api_solninfo(model.inner)[4]
    if stat == 1
        return MOI.FeasiblePoint
    else
        return MOI.UnknownResultStatus
    end
end

function LQOI.get_variable_primal_solution!(model::Optimizer, result)
    return c_api_getx(model.inner, result)
end

function LQOI.get_linear_primal_solution!(model::Optimizer, result)
    return c_api_getax(model.inner, result)
end

function LQOI.get_variable_dual_solution!(model::Optimizer, place)
    return c_api_getdj(model.inner, place)
end

function LQOI.get_linear_dual_solution!(model::Optimizer, place)
    return c_api_getpi(model.inner, place)
end

function LQOI.get_objective_value(model::Optimizer)
    return c_api_getobjval(model.inner)
end


text:
# Invertible CNN layer from Dinh et al. (2017)/Kingma & Dhariwal (2019)
# Author: Philipp Witte, pwitte3@gatech.edu
# Date: January 2020

using InvertibleNetworks, LinearAlgebra, Test

# Input
nx = 64
ny = 64
k = 20
n_in = 10
n_hidden = 20
batchsize = 2
type = "Flux"   # Flux CNN or residual block

# Input image
X = glorot_uniform(nx, ny, k, batchsize)
X0 = glorot_uniform(nx, ny, k, batchsize)

# 1x1 convolution and residual blocks
C = Conv1x1(k)

if type == "Flux"
    # Flux residual block (needs twice the number of output as input channels)
    model = Chain(
        Conv((3,3), n_in => n_hidden; pad=1),
        BatchNorm(n_hidden, relu),
        Conv((3,3), n_hidden => n_hidden; pad=1),
        BatchNorm(n_hidden, relu),
        Conv((3,3), n_hidden => n_hidden; pad=1),
        BatchNorm(n_hidden, relu)
    )
    RB = FluxBlock(model)
else
    RB = ResidualBlock(nx, ny, n_in, n_hidden, batchsize; fan=true)
end

# Invertible splitting layer
L = CouplingLayerGlow(C, RB; logdet=true)   # compute logdet

# Forward + backward
Y = L.forward(X)[1]
Y0, logdet = L.forward(X0)
ΔY = Y0 - Y
ΔX, X0_ = L.backward(ΔY, Y0)

@test isapprox(norm(X0_ - X0)/norm(X0), 0f0, atol=1f-2)


text:
# This file is part of Kpax3. License is MIT.

module TestPartition

  export Partition

  #=
  """
  # Partitions of n

  All the partitions of the set N = {1, ..., n} and their properties.

  ## Fields

  * `n` total number of units (integer to partition)
  * `B` total number of partitions of N (Bell number of `n`)
  * `C` total number of configuration classes
  * `k` total number of blocks within a particular configuration class
  * `partition` explicit partitions of the `n` units
  * `blocksize` block sizes corresponding to configuration classes
  * `blockcount` numbers of terms of various sizes
  * `index` indices of configuration classes within variable partition

  ## Details

  Let `k` represent the total number of blocks (clusters).

  Configuration classes are the number of distinct ways to partition the
  integer `n` into blocks of (n1, ..., nk) sizes, where

  * ng = #{i : i ∈ g},
  * n1 ≥ ... ng ≥ ... ≥ nk > 0
  * n1 + ... + nk = n.

  Partitions of `n` can also be represented by the numbers of terms of various
  sizes mj = #{g : ng = j} (j = 1, ..., n), where

  * m1 + ... mj + ... + mn = k
  * 1 * m1 + ... + j * mj + ... + n * mn = n
  """
  =#
  struct Partition
    n::Int
    B::Int
    C::Int
    k::Vector{Int}
    index::Vector{Int}
    blocksize::Matrix{Int}
    blockcount::Matrix{Int}
    partition::Matrix{Int}
  end

  #=
  """Partitions of 1"""
  =#
  function P1()
    n = 1
    B = 1
    C = 1

    k = ones(Int, C)
    index = ones(Int, C)

    blocksize  = ones(Int, n, C)
    blockcount = ones(Int, n, C)

    partition  = ones(Int, n, B)

    Partition(n, B, C, k, index, blocksize, blockcount, partition)
  end

  #=
  """Partitions of 2"""
  =#
  function P2()
    n = 2
    B = 2
    C = 2

    k = zeros(Int, C)
    index = zeros(Int, C)

    blocksize  = zeros(Int, n, C)
    blockcount = zeros(Int, n, C)

    partition  = zeros(Int, n, B)

    i = 0
    j = 0

    # k = 1

    # **
    j += 1
    index[j] = i + 1
    k[j] = 1
    blocksize[:,  j] .= [2, 0]
    blockcount[:, j] .= [0, 1]

    partition[:, i += 1] .= [1, 1]

    # k = 2

    # *|*
    j += 1
    index[j] = i + 1
    k[j] = 2
    blocksize[:,  j] .= [1, 1]
    blockcount[:, j] .= [2, 0]

    partition[:, i += 1] .= [1, 2]

    Partition(n, B, C, k, index, blocksize, blockcount, partition)
  end

  #=
  """Partitions of 3"""
  =#
  function P3()
    n = 3
    B = 5
    C = 3

    k = zeros(Int, C)
    index = zeros(Int, C)

    blocksize  = zeros(Int, n, C)
    blockcount = zeros(Int, n, C)

    partition  = zeros(Int, n, B)

    i = 0
    j = 0

    # k = 1

    # ***
    j += 1
    index[j] = i + 1
    k[j] = 1
    blocksize[:,  j] .= [3, 0, 0]
    blockcount[:, j] .= [0, 0, 1]

    partition[:, i += 1] .= [1, 1, 1]

    # k = 2

    # **|*
    j += 1
    index[j] = i + 1
    k[j] = 2
    blocksize[:,  j] .= [2, 1, 0]
    blockcount[:, j] .= [1, 1, 0]

    partition[:, i += 1] .= [1, 1, 2]
    partition[:, i += 1] .= [1, 2, 1]
    partition[:, i += 1] .= [2, 1, 1]

    # k = 3

    # *|*|*
    j += 1
    index[j] = i + 1
    k[j] = 3
    blocksize[:,  j] .= [1, 1, 1]
    blockcount[:, j] .= [3, 0, 0]

    partition[:, i += 1] .= [1, 2, 3]

    Partition(n, B, C, k, index, blocksize, blockcount, partition)
  end

  #=
  """Partitions of 4"""
  =#
  function P4()
    n = 4
    B = 15
    C = 5

    k = zeros(Int, C)
    index = zeros(Int, C)

    blocksize  = zeros(Int, n, C)
    blockcount = zeros(Int, n, C)

    partition  = zeros(Int, n, B)

    i = 0
    j = 0

    # k = 1

    # ****
    j += 1
    index[j] = i + 1
    k[j] = 1
    blocksize[:,  j] .= [4, 0, 0, 0]
    blockcount[:, j] .= [0, 0, 0, 1]

    partition[:, i += 1] .= [1, 1, 1, 1]

    # k = 2

    # ***|*
    j += 1
    index[j] = i + 1
    k[j] = 2
    blocksize[:,  j] .= [3, 1, 0, 0]
    blockcount[:, j] .= [1, 0, 1, 0]

    partition[:, i += 1] .= [1, 1, 1, 2]
    partition[:, i += 1] .= [1, 1, 2, 1]
    partition[:, i += 1] .= [1, 2, 1, 1]
    partition[:, i += 1] .= [2, 1, 1, 1]

    # **|**
    j += 1
    index[j] = i + 1
    k[j] = 2
    blocksize[:,  j] .= [2, 2, 0, 0]
    blockcount[:, j] .= [0, 2, 0, 0]

    partition[:, i += 1] .= [1, 1, 2, 2]
    partition[:, i += 1] .= [1, 2, 1, 2]
    partition[:, i += 1] .= [1, 2, 2, 1]

    # k = 3

    # **|*|*
    j += 1
    index[j] = i + 1
    k[j] = 3
    blocksize[:,  j] .= [2, 1, 1, 0]
    blockcount[:, j] .= [2, 1, 0, 0]

    partition[:, i += 1] .= [1, 1, 2, 3]
    partition[:, i += 1] .= [1, 2, 1, 3]
    partition[:, i += 1] .= [1, 2, 3, 1]
    partition[:, i += 1] .= [2, 1, 1, 3]
    partition[:, i += 1] .= [2, 1, 3, 1]
    partition[:, i += 1] .= [2, 3, 1, 1]

    # k = 4

    # *|*|*|*
    j += 1
    index[j] = i + 1
    k[j] = 4
    blocksize[:,  j] .= [1, 1, 1, 1]
    blockcount[:, j] .= [4, 0, 0, 0]

    partition[:, i += 1] .= [1, 2, 3, 4]

    Partition(n, B, C, k, index, blocksize, blockcount, partition)
  end

  #=
  """Partitions of 5"""
  =#
  function P5()
    n = 5
    B = 52
    C = 7

    k = zeros(Int, C)
    index = zeros(Int, C)

    blocksize  = zeros(Int, n, C)
    blockcount = zeros(Int, n, C)

    partition  = zeros(Int, n, B)

    i = 0
    j = 0

    # k = 1

    # *****
    j += 1
    index[j] = i + 1
    k[j] = 1
    blocksize[:,  j] .= [5, 0, 0, 0, 0]
    blockcount[:, j] .= [0, 0, 0, 0, 1]

    partition[:, i += 1] .= [1, 1, 1, 1, 1]

    # k = 2

    # ****|*
    j += 1
    index[j] = i + 1
    k[j] = 2
    blocksize[:,  j] .= [4, 1, 0, 0, 0]
    blockcount[:, j] .= [1, 0, 0, 1, 0]

    partition[:, i += 1] .= [1, 1, 1, 1, 2]
    partition[:, i += 1] .= [1, 1, 1, 2, 1]
    partition[:, i += 1] .= [1, 1, 2, 1, 1]
    partition[:, i += 1] .= [1, 2, 1, 1, 1]
    partition[:, i += 1] .= [2, 1, 1, 1, 1]

    # ***|**
    j += 1
    index[j] = i + 1
    k[j] = 2
    blocksize[:,  j] .= [3, 2, 0, 0, 0]
    blockcount[:, j] .= [0, 1, 1, 0, 0]

    partition[:, i += 1] .= [1, 1, 1, 2, 2]
    partition[:, i += 1] .= [1, 1, 2, 1, 2]
    partition[:, i += 1] .= [1, 2, 1, 1, 2]
    partition[:, i += 1] .= [2, 1, 1, 1, 2]
    partition[:, i += 1] .= [1, 1, 2, 2, 1]
    partition[:, i += 1] .= [1, 2, 1, 2, 1]
    partition[:, i += 1] .= [2, 1, 1, 2, 1]
    partition[:, i += 1] .= [1, 2, 2, 1, 1]
    partition[:, i += 1] .= [2, 1, 2, 1, 1]
    partition[:, i += 1] .= [2, 2, 1, 1, 1]

    # k = 3

    # **|**|*
    j += 1
    index[j] = i + 1
    k[j] = 3
    blocksize[:,  j] .= [2, 2, 1, 0, 0]
    blockcount[:, j] .= [1, 2, 0, 0, 0]

    partition[:, i += 1] .= [1, 1, 2, 2, 3]
    partition[:, i += 1] .= [1, 1, 2, 3, 2]
    partition[:, i += 1] .= [1, 1, 3, 2, 2]
    partition[:, i += 1] .= [1, 2, 1, 2, 3]
    partition[:, i += 1] .= [1, 2, 1, 3, 2]
    partition[:, i += 1] .= [1, 3, 1, 2, 2]
    partition[:, i += 1] .= [1, 2, 2, 1, 3]
    partition[:, i += 1] .= [1, 2, 3, 1, 2]
    partition[:, i += 1] .= [1, 3, 2, 1, 2]
    partition[:, i += 1] .= [1, 2, 2, 3, 1]
    partition[:, i += 1] .= [1, 2, 3, 2, 1]
    partition[:, i += 1] .= [1, 3, 2, 2, 1]
    partition[:, i += 1] .= [3, 1, 1, 2, 2]
    partition[:, i += 1] .= [3, 1, 2, 1, 2]
    partition[:, i += 1] .= [3, 1, 2, 2, 1]

    # ***|*|*
    j += 1
    index[j] = i + 1
    k[j] = 3
    blocksize[:,  j] .= [3, 1, 1, 0, 0]
    blockcount[:, j] .= [2, 0, 1, 0, 0]

    partition[:, i += 1] .= [1, 1, 1, 2, 3]
    partition[:, i += 1] .= [1, 1, 2, 1, 3]
    partition[:, i += 1] .= [1, 1, 2, 3, 1]
    partition[:, i += 1] .= [1, 2, 1, 1, 3]
    partition[:, i += 1] .= [1, 2, 1, 3, 1]
    partition[:, i += 1] .= [1, 2, 3, 1, 1]
    partition[:, i += 1] .= [2, 1, 1, 1, 3]
    partition[:, i += 1] .= [2, 1, 1, 3, 1]
    partition[:, i += 1] .= [2, 1, 3, 1, 1]
    partition[:, i += 1] .= [2, 3, 1, 1, 1]

    # k = 4

    # **|*|*|*
    j += 1
    index[j] = i + 1
    k[j] = 4
    blocksize[:,  j] .= [2, 1, 1, 1, 0]
    blockcount[:, j] .= [3, 1, 0, 0, 0]

    partition[:, i += 1] .= [1, 1, 2, 3, 4]
    partition[:, i += 1] .= [1, 2, 1, 3, 4]
    partition[:, i += 1] .= [1, 2, 3, 1, 4]
    partition[:, i += 1] .= [1, 2, 3, 4, 1]
    partition[:, i += 1] .= [2, 1, 1, 3, 4]
    partition[:, i += 1] .= [2, 1, 3, 1, 4]
    partition[:, i += 1] .= [2, 1, 3, 4, 1]
    partition[:, i += 1] .= [2, 3, 1, 1, 4]
    partition[:, i += 1] .= [2, 3, 1, 4, 1]
    partition[:, i += 1] .= [2, 3, 4, 1, 1]

    # k = 5

    # *|*|*|*|*
    j += 1
    index[j] = i + 1
    k[j] = 5
    blocksize[:,  j] .= [1, 1, 1, 1, 1]
    blockcount[:, j] .= [5, 0, 0, 0, 0]

    partition[:, i += 1] .= [1, 2, 3, 4, 5]

    Partition(n, B, C, k, index, blocksize, blockcount, partition)
  end

  #=
  """Partitions of 6"""
  =#
  function P6()
    n = 6
    B = 203
    C = 11

    k = zeros(Int, C)
    index = zeros(Int, C)

    blocksize  = zeros(Int, n, C)
    blockcount = zeros(Int, n, C)

    partition  = zeros(Int, n, B)

    i = 0
    j = 0

    # k = 1

    # ******
    j += 1
    index[j] = i + 1
    k[j] = 1
    blocksize[:,  j] .= [6, 0, 0, 0, 0, 0]
    blockcount[:, j] .= [0, 0, 0, 0, 0, 6]

    partition[:, i += 1] .= [1, 1, 1, 1, 1, 1]

    # k = 2

    # *****|*
    j += 1
    index[j] = i + 1
    k[j] = 2
    blocksize[:,  j] .= [5, 1, 0, 0, 0, 0]
    blockcount[:, j] .= [1, 0, 0, 0, 1, 0]

    partition[:, i += 1] .= [1, 1, 1, 1, 1, 2]
    partition[:, i += 1] .= [1, 1, 1, 1, 2, 1]
    partition[:, i += 1] .= [1, 1, 1, 2, 1, 1]
    partition[:, i += 1] .= [1, 1, 2, 1, 1, 1]
    partition[:, i += 1] .= [1, 2, 1, 1, 1, 1]
    partition[:, i += 1] .= [2, 1, 1, 1, 1, 1]

    # ****|**
    j += 1
    index[j] = i + 1
    k[j] = 2
    blocksize[:,  j] .= [4, 2, 0, 0, 0, 0]
    blockcount[:, j] .= [0, 1, 0, 1, 0, 0]

    partition[:, i += 1] .= [1, 1, 1, 1, 2, 2]
    partition[:, i += 1] .= [1, 1, 1, 2, 1, 2]
    partition[:, i += 1] .= [1, 1, 1, 2, 2, 1]
    partition[:, i += 1] .= [1, 1, 2, 1, 1, 2]
    partition[:, i += 1] .= [1, 1, 2, 1, 2, 1]
    partition[:, i += 1] .= [1, 1, 2, 2, 1, 1]
    partition[:, i += 1] .= [1, 2, 1, 1, 1, 2]
    partition[:, i += 1] .= [1, 2, 1, 1, 2, 1]
    partition[:, i += 1] .= [1, 2, 1, 2, 1, 1]
    partition[:, i += 1] .= [1, 2, 2, 1, 1, 1]
    partition[:, i += 1] .= [2, 1, 1, 1, 1, 2]
    partition[:, i += 1] .= [2, 1, 1, 1, 2, 1]
    partition[:, i += 1] .= [2, 1, 1, 2, 1, 1]
    partition[:, i += 1] .= [2, 1, 2, 1, 1, 1]
    partition[:, i += 1] .= [2, 2, 1, 1, 1, 1]

    # ***|***
    j += 1
    index[j] = i + 1
    k[j] = 2
    blocksize[:,  j] .= [3, 3, 0, 0, 0, 0]
    blockcount[:, j] .= [0, 0, 2, 0, 0, 0]

    partition[:, i += 1] .= [1, 1, 1, 2, 2, 2]
    partition[:, i += 1] .= [1, 1, 2, 1, 2, 2]
    partition[:, i += 1] .= [1, 1, 2, 2, 1, 2]
    partition[:, i += 1] .= [1, 1, 2, 2, 2, 1]
    partition[:, i += 1] .= [1, 2, 1, 1, 2, 2]
    partition[:, i += 1] .= [1, 2, 1, 2, 1, 2]
    partition[:, i += 1] .= [1, 2, 1, 2, 2, 1]
    partition[:, i += 1] .= [1, 2, 2, 1, 1, 2]
    partition[:, i += 1] .= [1, 2, 2, 1, 2, 1]
    partition[:, i += 1] .= [1, 2, 2, 2, 1, 1]

    # k = 3

    # ****|*|*
    j += 1
    index[j] = i + 1
    k[j] = 3
    blocksize[:,  j] .= [4, 1, 1, 0, 0, 0]
    blockcount[:, j] .= [2, 0, 0, 1, 0, 0]

    partition[:, i += 1] .= [1, 1, 1, 1, 2, 3]
    partition[:, i += 1] .= [1, 1, 1, 2, 1, 3]
    partition[:, i += 1] .= [1, 1, 1, 2, 3, 1]
    partition[:, i += 1] .= [1, 1, 2, 1, 1, 3]
    partition[:, i += 1] .= [1, 1, 2, 1, 3, 1]
    partition[:, i += 1] .= [1, 1, 2, 3, 1, 1]
    partition[:, i += 1] .= [1, 2, 1, 1, 1, 3]
    partition[:, i += 1] .= [1, 2, 1, 1, 3, 1]
    partition[:, i += 1] .= [1, 2, 1, 3, 1, 1]
    partition[:, i += 1] .= [1, 2, 3, 1, 1, 1]
    partition[:, i += 1] .= [2, 1, 1, 1, 1, 3]
    partition[:, i += 1] .= [2, 1, 1, 1, 3, 1]
    partition[:, i += 1] .= [2, 1, 1, 3, 1, 1]
    partition[:, i += 1] .= [2, 1, 3, 1, 1, 1]
    partition[:, i += 1] .= [2, 3, 1, 1, 1, 1]

    # ***|**|*
    j += 1
    index[j] = i + 1
    k[j] = 3
    blocksize[:,  j] .= [3, 2, 1, 0, 0, 0]
    blockcount[:, j] .= [1, 1, 1, 0, 0, 0]

    partition[:, i += 1] .= [1, 1, 1, 2, 2, 3]
    partition[:, i += 1] .= [1, 1, 1, 2, 3, 2]
    partition[:, i += 1] .= [1, 1, 1, 3, 2, 2]
    partition[:, i += 1] .= [1, 1, 2, 1, 2, 3]
    partition[:, i += 1] .= [1, 1, 2, 1, 3, 2]
    partition[:, i += 1] .= [1, 1, 3, 1, 2, 2]
    partition[:, i += 1] .= [1, 1, 2, 2, 1, 3]
    partition[:, i += 1] .= [1, 1, 2, 3, 1, 2]
    partition[:, i += 1] .= [1, 1, 3, 2, 1, 2]
    partition[:, i += 1] .= [1, 1, 2, 2, 3, 1]
    partition[:, i += 1] .= [1, 1, 2, 3, 2, 1]
    partition[:, i += 1] .= [1, 1, 3, 2, 2, 1]
    partition[:, i += 1] .= [1, 2, 1, 1, 2, 3]
    partition[:, i += 1] .= [1, 2, 1, 1, 3, 2]
    partition[:, i += 1] .= [1, 3, 1, 1, 2, 2]
    partition[:, i += 1] .= [1, 2, 1, 2, 1, 3]
    partition[:, i += 1] .= [1, 2, 1, 3, 1, 2]
    partition[:, i += 1] .= [1, 3, 1, 2, 1, 2]
    partition[:, i += 1] .= [1, 2, 1, 2, 3, 1]
    partition[:, i += 1] .= [1, 2, 1, 3, 2, 1]
    partition[:, i += 1] .= [1, 3, 1, 2, 2, 1]
    partition[:, i += 1] .= [1, 2, 2, 1, 1, 3]
    partition[:, i += 1] .= [1, 2, 3, 1, 1, 2]
    partition[:, i += 1] .= [1, 3, 2, 1, 1, 2]
    partition[:, i += 1] .= [1, 2, 2, 1, 3, 1]
    partition[:, i += 1] .= [1, 2, 3, 1, 2, 1]
    partition[:, i += 1] .= [1, 3, 2, 1, 2, 1]
    partition[:, i += 1] .= [1, 2, 2, 3, 1, 1]
    partition[:, i += 1] .= [1, 2, 3, 2, 1, 1]
    partition[:, i += 1] .= [1, 3, 2, 2, 1, 1]
    partition[:, i += 1] .= [2, 1, 1, 1, 2, 3]
    partition[:, i += 1] .= [2, 1, 1, 1, 3, 2]
    partition[:, i += 1] .= [3, 1, 1, 1, 2, 2]
    partition[:, i += 1] .= [2, 1, 1, 2, 1, 3]
    partition[:, i += 1] .= [2, 1, 1, 3, 1, 2]
    partition[:, i += 1] .= [3, 1, 1, 2, 1, 2]
    partition[:, i += 1] .= [2, 1, 1, 2, 3, 1]
    partition[:, i += 1] .= [2, 1, 1, 3, 2, 1]
    partition[:, i += 1] .= [3, 1, 1, 2, 2, 1]
    partition[:, i += 1] .= [2, 1, 2, 1, 1, 3]
    partition[:, i += 1] .= [2, 1, 3, 1, 1, 2]
    partition[:, i += 1] .= [3, 1, 2, 1, 1, 2]
    partition[:, i += 1] .= [2, 1, 2, 1, 3, 1]
    partition[:, i += 1] .= [2, 1, 3, 1, 2, 1]
    partition[:, i += 1] .= [3, 1, 2, 1, 2, 1]
    partition[:, i += 1] .= [2, 1, 2, 3, 1, 1]
    partition[:, i += 1] .= [2, 1, 3, 2, 1, 1]
    partition[:, i += 1] .= [3, 1, 2, 2, 1, 1]
    partition[:, i += 1] .= [2, 2, 1, 1, 1, 3]
    partition[:, i += 1] .= [2, 3, 1, 1, 1, 2]
    partition[:, i += 1] .= [3, 2, 1, 1, 1, 2]
    partition[:, i += 1] .= [2, 2, 1, 1, 3, 1]
    partition[:, i += 1] .= [2, 3, 1, 1, 2, 1]
    partition[:, i += 1] .= [3, 2, 1, 1, 2, 1]
    partition[:, i += 1] .= [2, 2, 1, 3, 1, 1]
    partition[:, i += 1] .= [2, 3, 1, 2, 1, 1]
    partition[:, i += 1] .= [3, 2, 1, 2, 1, 1]
    partition[:, i += 1] .= [2, 2, 3, 1, 1, 1]
    partition[:, i += 1] .= [2, 3, 2, 1, 1, 1]
    partition[:, i += 1] .= [3, 2, 2, 1, 1, 1]

    # **|**|**
    j += 1
    index[j] = i + 1
    k[j] = 3
    blocksize[:,  j] .= [2, 2, 2, 0, 0, 0]
    blockcount[:, j] .= [0, 3, 0, 0, 0, 0]

    partition[:, i += 1] .= [1, 1, 2, 2, 3, 3]
    partition[:, i += 1] .= [1, 1, 2, 3, 2, 3]
    partition[:, i += 1] .= [1, 1, 2, 3, 3, 2]
    partition[:, i += 1] .= [1, 2, 1, 2, 3, 3]
    partition[:, i += 1] .= [1, 2, 1, 3, 2, 3]
    partition[:, i += 1] .= [1, 2, 1, 3, 3, 2]
    partition[:, i += 1] .= [1, 2, 2, 1, 3, 3]
    partition[:, i += 1] .= [1, 2, 3, 1, 2, 3]
    partition[:, i += 1] .= [1, 2, 3, 1, 3, 2]
    partition[:, i += 1] .= [1, 2, 2, 3, 1, 3]
    partition[:, i += 1] .= [1, 2, 3, 2, 1, 3]
    partition[:, i += 1] .= [1, 2, 3, 3, 1, 2]
    partition[:, i += 1] .= [1, 2, 2, 3, 3, 1]
    partition[:, i += 1] .= [1, 2, 3, 2, 3, 1]
    partition[:, i += 1] .= [1, 2, 3, 3, 2, 1]

    # k = 4

    # ***|*|*|*
    j += 1
    index[j] = i + 1
    k[j] = 4
    blocksize[:,  j] .= [3, 1, 1, 1, 0, 0]
    blockcount[:, j] .= [3, 0, 1, 0, 0, 0]

    partition[:, i += 1] .= [1, 1, 1, 2, 3, 4]
    partition[:, i += 1] .= [1, 1, 2, 1, 3, 4]
    partition[:, i += 1] .= [1, 1, 2, 3, 1, 4]
    partition[:, i += 1] .= [1, 1, 2, 3, 4, 1]
    partition[:, i += 1] .= [1, 2, 1, 1, 3, 4]
    partition[:, i += 1] .= [1, 2, 1, 3, 1, 4]
    partition[:, i += 1] .= [1, 2, 1, 3, 4, 1]
    partition[:, i += 1] .= [1, 2, 3, 1, 1, 4]
    partition[:, i += 1] .= [1, 2, 3, 1, 4, 1]
    partition[:, i += 1] .= [1, 2, 3, 4, 1, 1]
    partition[:, i += 1] .= [2, 1, 1, 1, 3, 4]
    partition[:, i += 1] .= [2, 1, 1, 3, 1, 4]
    partition[:, i += 1] .= [2, 1, 1, 3, 4, 1]
    partition[:, i += 1] .= [2, 1, 3, 1, 1, 4]
    partition[:, i += 1] .= [2, 1, 3, 1, 4, 1]
    partition[:, i += 1] .= [2, 1, 3, 4, 1, 1]
    partition[:, i += 1] .= [2, 3, 1, 1, 1, 4]
    partition[:, i += 1] .= [2, 3, 1, 1, 4, 1]
    partition[:, i += 1] .= [2, 3, 1, 4, 1, 1]
    partition[:, i += 1] .= [2, 3, 4, 1, 1, 1]

    # **|**|*|*
    j += 1
    index[j] = i + 1
    k[j] = 4
    blocksize[:,  j] .= [2, 2, 1, 1, 0, 0]
    blockcount[:, j] .= [2, 2, 0, 0, 0, 0]

    partition[:, i += 1] .= [1, 1, 2, 2, 3, 4]
    partition[:, i += 1] .= [1, 1, 2, 3, 2, 4]
    partition[:, i += 1] .= [1, 1, 2, 3, 4, 2]
    partition[:, i += 1] .= [1, 1, 3, 2, 2, 4]
    partition[:, i += 1] .= [1, 1, 3, 2, 4, 2]
    partition[:, i += 1] .= [1, 1, 3, 4, 2, 2]
    partition[:, i += 1] .= [1, 2, 1, 2, 3, 4]
    partition[:, i += 1] .= [1, 2, 1, 3, 2, 4]
    partition[:, i += 1] .= [1, 2, 1, 3, 4, 2]
    partition[:, i += 1] .= [1, 3, 1, 2, 2, 4]
    partition[:, i += 1] .= [1, 3, 1, 2, 4, 2]
    partition[:, i += 1] .= [1, 3, 1, 4, 2, 2]
    partition[:, i += 1] .= [1, 2, 2, 1, 3, 4]
    partition[:, i += 1] .= [1, 2, 3, 1, 2, 4]
    partition[:, i += 1] .= [1, 2, 3, 1, 4, 2]
    partition[:, i += 1] .= [1, 3, 2, 1, 2, 4]
    partition[:, i += 1] .= [1, 3, 2, 1, 4, 2]
    partition[:, i += 1] .= [1, 3, 4, 1, 2, 2]
    partition[:, i += 1] .= [1, 2, 2, 3, 1, 4]
    partition[:, i += 1] .= [1, 2, 3, 2, 1, 4]
    partition[:, i += 1] .= [1, 2, 3, 4, 1, 2]
    partition[:, i += 1] .= [1, 3, 2, 2, 1, 4]
    partition[:, i += 1] .= [1, 3, 2, 4, 1, 2]
    partition[:, i += 1] .= [1, 3, 4, 2, 1, 2]
    partition[:, i += 1] .= [1, 2, 2, 3, 4, 1]
    partition[:, i += 1] .= [1, 2, 3, 2, 4, 1]
    partition[:, i += 1] .= [1, 2, 3, 4, 2, 1]
    partition[:, i += 1] .= [1, 3, 2, 2, 4, 1]
    partition[:, i += 1] .= [1, 3, 2, 4, 2, 1]
    partition[:, i += 1] .= [1, 3, 4, 2, 2, 1]
    partition[:, i += 1] .= [3, 1, 1, 2, 2, 4]
    partition[:, i += 1] .= [3, 1, 1, 2, 4, 2]
    partition[:, i += 1] .= [3, 1, 1, 4, 2, 2]
    partition[:, i += 1] .= [3, 1, 2, 1, 2, 4]
    partition[:, i += 1] .= [3, 1, 2, 1, 4, 2]
    partition[:, i += 1] .= [3, 1, 4, 1, 2, 2]
    partition[:, i += 1] .= [3, 1, 2, 2, 1, 4]
    partition[:, i += 1] .= [3, 1, 2, 4, 1, 2]
    partition[:, i += 1] .= [3, 1, 4, 2, 1, 2]
    partition[:, i += 1] .= [3, 1, 2, 2, 4, 1]
    partition[:, i += 1] .= [3, 1, 2, 4, 2, 1]
    partition[:, i += 1] .= [3, 1, 4, 2, 2, 1]
    partition[:, i += 1] .= [3, 4, 1, 1, 2, 2]
    partition[:, i += 1] .= [3, 4, 1, 2, 1, 2]
    partition[:, i += 1] .= [3, 4, 1, 2, 2, 1]

    # k = 5

    # **|*|*|*|*
    j += 1
    index[j] = i + 1
    k[j] = 5
    blocksize[:,  j] .= [2, 1, 1, 1, 1, 0]
    blockcount[:, j] .= [4, 1, 0, 0, 0, 0]

    partition[:, i += 1] .= [1, 1, 2, 3, 4, 5]
    partition[:, i += 1] .= [1, 2, 1, 3, 4, 5]
    partition[:, i += 1] .= [1, 2, 3, 1, 4, 5]
    partition[:, i += 1] .= [1, 2, 3, 4, 1, 5]
    partition[:, i += 1] .= [1, 2, 3, 4, 5, 1]
    partition[:, i += 1] .= [2, 1, 1, 3, 4, 5]
    partition[:, i += 1] .= [2, 1, 3, 1, 4, 5]
    partition[:, i += 1] .= [2, 1, 3, 4, 1, 5]
    partition[:, i += 1] .= [2, 1, 3, 4, 5, 1]
    partition[:, i += 1] .= [2, 3, 1, 1, 4, 5]
    partition[:, i += 1] .= [2, 3, 1, 4, 1, 5]
    partition[:, i += 1] .= [2, 3, 1, 4, 5, 1]
    partition[:, i += 1] .= [2, 3, 4, 1, 1, 5]
    partition[:, i += 1] .= [2, 3, 4, 1, 5, 1]
    partition[:, i += 1] .= [2, 3, 4, 5, 1, 1]

    # k = 6

    # *|*|*|*|*|*
    j += 1
    index[j] = i + 1
    k[j] = 6
    blocksize[:,  j] .= [1, 1, 1, 1, 1, 1]
    blockcount[:, j] .= [6, 0, 0, 0, 0, 0]

    partition[:, i += 1] .= [1, 2, 3, 4, 5, 6]

    Partition(n, B, C, k, index, blocksize, blockcount, partition)
  end

  function Partition(n::Int)
    if n == 1
      P1()
    elseif n == 2
      P2()
    elseif n == 3
      P3()
    elseif n == 4
      P4()
    elseif n == 5
      P5()
    elseif n == 6
      P6()
    else
      throw(KDomainError("Argument n must be positive and lesser than 7."))
    end
  end
end


text:
# SCHEME
abstract type Scheme end

# By default, scheme definition can ignore time step
(scheme::Scheme)(flux, mesh, w, i_face, dt) = scheme(flux, mesh, w, i_face)

####################################
struct Centered <: Scheme end
####################################

function (::Centered)(flux::LinearAdvectionFlux, mesh, w, i_face)
    n = normal_vector(mesh, i_face)
    i_cell_1, i_cell_2 = cells_next_to_inner_face(mesh, i_face)
    return 0.5*(flux(w[i_cell_1], n) + flux(w[i_cell_2], n))
end

####################################
struct Downwind <: Scheme end
####################################

function (::Downwind)(flux::LinearAdvectionFlux, mesh, w, i_face)
    n = normal_vector(mesh, i_face)
    i_cell_1, i_cell_2 = cells_next_to_inner_face(mesh, i_face)
    λ = eigvals(flux, (w[i_cell_1] + w[i_cell_2])/2, n)
    return λ < zero(λ) ? flux(w[i_cell_1], n) : flux(w[i_cell_2], n)
end

##################################
struct Upwind <: Scheme end
##################################

# Single-wave problems
function (::Upwind)(flux::LinearAdvectionFlux, mesh, w, i_face)
    n = normal_vector(mesh, i_face)
    i_cell_1, i_cell_2 = cells_next_to_inner_face(mesh, i_face)
    λ = eigvals(flux, (w[i_cell_1] + w[i_cell_2])/2, n)
    return λ > zero(λ) ? flux(w[i_cell_1], n) : flux(w[i_cell_2], n)
end

# Multi-wave problems
function (::Upwind)(flux, mesh, w, i_face)
    n = normal_vector(mesh, i_face)
    i_cell_1, i_cell_2 = cells_next_to_inner_face(mesh, i_face)
    w_mean = (w[i_cell_1] + w[i_cell_2])/2
    λ, R = eigen(flux, w_mean, n)
    left_to_right = λ .> zero(λ)
    if all(left_to_right)
        return flux(w[i_cell_1], n)
    elseif !any(left_to_right)
        return flux(w[i_cell_2], n)
    else
        f₁ = flux(w[i_cell_1], n)
        f₂ = flux(w[i_cell_2], n)
        decentrement = R \ (f₁ - f₂)
        return (f₁ + f₂ + R*diagm(sign.(λ))*decentrement)/2
        # L = inv(R)
        # return (f₁ + f₂ + R*diagm(sign.(λ))*L*(f₁ - f₂))/2
    end
end

###################################
abstract type BoundaryCondition end
###################################

(scheme::BoundaryCondition)(model, mesh, w, i_face, dt) = scheme(model, mesh, w, i_face)

#########################################
struct NeumannBC <: BoundaryCondition end
#########################################

function (::NeumannBC)(flux, mesh, w, i_face)
    i_cell_1 = cell_next_to_boundary_face(mesh, i_face)
    return flux(w[i_cell_1], normal_vector(mesh, i_face))
end

################################################################################
#                                     Div                                      #
################################################################################

numerical_flux(flux, mesh, w, Φ, i_face, dt) = Φ(flux, mesh, w, i_face, dt)

function div!(Δw, flux, mesh, w, scheme::Scheme, dt=0.0)
    @inbounds for i_face in inner_faces(mesh)
        ϕ = numerical_flux(flux, mesh, w, scheme, i_face, dt)
        i_cell_1, i_cell_2 = cells_next_to_inner_face(mesh, i_face)
        Δw[i_cell_1] += ϕ * face_area(mesh, i_face) / cell_volume(mesh, i_cell_1)
        Δw[i_cell_2] -= ϕ * face_area(mesh, i_face) / cell_volume(mesh, i_cell_2)
    end
end

function div!(Δw, flux, mesh, w, scheme::BoundaryCondition, dt=0.0)
    @inbounds for i_face in boundary_faces(mesh)
        ϕ = numerical_flux(flux, mesh, w, scheme, i_face, dt)
        i_cell = cell_next_to_boundary_face(mesh, i_face)
        Δw[i_cell] += ϕ * face_area(mesh, i_face) / cell_volume(mesh, i_cell)
    end
end

function div!(Δw, flux, mesh, w, schemes::NTuple{N, Union{Scheme, BoundaryCondition}}, dt=0.0) where N
    for scheme in schemes
        div!(Δw, flux, mesh, w, scheme, dt)
    end
end

flux_in_cell(flux, mesh, w, scheme, dt, i_face, i_cell) = numerical_flux(flux, mesh, w, scheme, i_face, dt) * face_area(mesh, i_face) / cell_volume(mesh, i_cell)
function Δw_type(flux, mesh, w, schemes::NTuple{N, Union{Scheme, BoundaryCondition}}, dt) where N
    T = Base.return_types(flux_in_cell, typeof.((flux, mesh, w, schemes[1], dt, first(inner_faces(mesh)), first(all_cells(mesh)))))[1]
    if T in (Union{}, Any)
        return eltype(w)
    else
        return T
    end
end
Δw_type(f, m, w, s::Union{Scheme, BoundaryCondition}, dt) = Δw_type(f, m, w, (s,), dt)

function div(flux, mesh, w, schemes::NTuple{N, Union{Scheme, BoundaryCondition}}, dt=0.0) where N
    Δw = zeros(Δw_type(flux, mesh, w, schemes, dt), size(w))
    div!(Δw, flux, mesh, w, schemes, dt)
    return Δw
end

div(flux, mesh, w, scheme::Union{Scheme, BoundaryCondition}, dt=0.0) = div(flux, mesh, w, (scheme,), dt)