text:
## deriv.jl - derivative of a single ExNode

const DIFF_PHS = Set([:w, :x, :y, :z, :i, :j, :k,])

isparameters(a) = isa(a, Expr) && a.head == :parameters

function without_types(pat)
    rpat = copy(pat)
    for i=2:length(rpat.args)
        a = rpat.args[i]
        if !isparameters(a)  # parameters aren't modified
            rpat.args[i] = isa(a, Expr) ? a.args[1] : a
        end
    end
    return rpat
end

function get_arg_names(pat)
    return [isa(a, Expr) ? a.args[1] : a for a in pat.args[2:end] if !isparameters(a)]
end

function get_arg_types(pat)
    return [isa(a, Expr) ? Core.eval(Base, a.args[2]) : Any
            for a in pat.args[2:end] if !isparameters(a)]
end

function match_rule(rule, ex, dep_vals, idx)
    tpat, (vname, rpat) = rule
    vidx = findfirst(isequal(vname), get_arg_names(tpat))
    if idx != vidx
        return nothing
    end
    dep_types = get_arg_types(tpat)
    if length(dep_vals) != length(dep_types) ||
        !all(isa(v, t) for (v, t) in zip(dep_vals, dep_types))
        return nothing
    end
    pat = without_types(tpat)
    ex_ = without_keywords(ex)
    if !matchingex(pat, ex_; phs=DIFF_PHS)
        return nothing
    else
        return pat => rpat
    end
end

function rewrite_with_keywords(ex, pat, rpat)
    if rpat isa Expr && rpat.head == :call
        op, args, kw_args = parse_call_expr(ex)
        ex_no_kw = make_call_expr(op, args)
        rex_no_kw = rewrite(ex_no_kw, pat, rpat; phs=DIFF_PHS)
        rex = with_keywords(rex_no_kw, kw_args)
    else
        rex = rewrite(ex, pat, rpat; phs=DIFF_PHS)
    end
    return rex
end

function deriv(ex, dep_vals, idx::Int)
    rex = nothing
    for rule in DIFF_RULES
        m = match_rule(rule, ex, dep_vals, idx)
        if m != nothing
            pat, rpat = m
            rex = rewrite_with_keywords(ex, pat, rpat)
            break
        end
    end
    if rex == nothing
        error("Can't find differentiation rule for $ex at $idx " *
              "with types $(map(typeof, dep_vals))")
    end
    return rex
end

"""
Internal function for finding rule by function name
"""
function find_rules_for(fun)
    return [r for r in DIFF_RULES if r[1].args[1] == fun]
end


text:
using MCMCDiagnostics
using Measurements
import StatsBase: mean_and_std

export expect

# Expected value of x log-weighted by ℓ
function expect(ℓ,x)
    ℓmax = maximum(ℓ)
    xw = 0.0
    ∑w = 0.0
    for (ℓi, xi) in zip(ℓ,x)
        wi = exp(ℓi - ℓmax)
        xw += xi * wi
        ∑w += wi
    end
    return xw / ∑w
end

function expect(x)
    (μ,σ) = mean_and_std(x)
    neff = effective_sample_size(x)
    return μ ± (σ / √neff)
end

export @cleanbreak

macro cleanbreak(ex)
    quote
        try
            $(esc(ex))
        catch e
            if e isa InterruptException
                @warn "Computation interrupted"
            else
                rethrow()
            end
        end
    end
end

ext: jl


text:
# To do: pass the storages f, ∆X, ∆F, Q, β.  Also, think about ways to reduce this memory.

# To use andersonaccel!, implement anderson_SALT! that accepts SALTSol as an initial guess
# and g! that takes SALTSol and returns SALTSol.  anderson_SALT! should create a version of g!
# that takes and returns vectors using CatViews.

# I will need to implement `reinterpret` for PETSc vectors to view complex PETSc vectors as
# a real PETSc vector.
function lsol2rvec(lsol::LasingSol)
    m_active = lsol.m_active
    ψr = reinterpret.(Ref(Float), @view(lsol.ψ[lsol.m_active]))
    # ψr = lsol.ψ[lsol.m_active]  # complex version

    # # For the Anderson in the complex domain
    # ωa² = lsol.ω[lsol.m_active] .+ im .* lsol.a²[lsol.m_active]
    # return CatView(ωa², lsol.ψ[lsol.m_active]...)

    # For the Anderson in the real domain
    return CatView(@view(lsol.ω[lsol.m_active]), @view(lsol.a²[lsol.m_active]), ψr...)

    # return CatView(ψr...)
    # return CatView(lsol.ω[lsol.m_active], lsol.a²[lsol.m_active])
end

# Create the fixed-point function to use in the Anderson acceleration.
function create_fpfun(lsol::LasingSol, lvar::LasingVar, gp::GainProfile)
    # When x is created from lsol by lsol2rvec, x points to the memory of lsol, so the update in
    # lsol is automatically reflected in x.
    function fpfun!(y, x)
        update_lsol!(lsol, lvar, gp)  # x is updated to g(xₖ) (but this x is not xₖ₊₁)

        # Uncomment the following two lines for the Anderson in the complex domain.
        # @info "This is complex Anderson."
        # y[1:length(lsol.m_active)] .= lsol.ω[lsol.m_active] .+ im .* lsol.a²[lsol.m_active]

        return nothing
    end

    return fpfun!
end

function create_normfun(lsol::LasingSol, lvar::LasingVar, gp::GainProfile, εc::AbsVecComplexF;
                        verbose::Bool=true, msgprefix::String="    ")
    function normfun!(x::AbsVec{<:Number}, k::Integer=0, lleq₀::Real=Inf)
        init_lvar!(lvar, lsol, gp, εc)
        lleq = norm_leq(lsol, lvar, gp)
        if verbose
            k==0 ? println(msgprefix * "Initial residual norm: ‖leq₀‖ = $lleq") :
                   println(msgprefix * "k = $k: ‖leq‖/‖leq₀‖ = $(lleq/lleq₀)")
        end

        return lleq
    end

    return normfun!
end

# Solve the fixed-point equation g(x) = x to given relative and absolute tolerances,
# returning the solution x, via Anderson acceleration of a fixed-point iteration starting at
# a given initial x (which should be of the correct type to hold the result).
#
# x should be set x = lsol2rvec(lsol) outside.
#
# This implementation operates in-place as much as possible.
# This function will be removed if an external implementation of the Anderson acceleration is used.
function anderson!(g!::Any,  # called by g!(y,x); overwrites a given output vector y with g(x)
                   x::AbsVec{T};  # initial guess; overwritten at each iteration step (note g! itself does not overwrite x)
                   norm=(x::AbsVec{T}, k::Integer, l₀::Real)->norm(x), # stopping criterion: norm(Δx) ≤ atol + rtol*norm(x)
                   m::Integer=2, # number of additional x's kept in algorithm; m = 0 means unaccelerated iteration
                   rtol::Real=τᵣ₀,  # relative tolerance; consider using τᵣ₀
                   atol::Real=τₐ₀,  # absolute tolerance
                   maxit::Int=typemax(Int)) where {T<:Number}  # number of maximum number of iteration steps
    m ≥ 0 || @error "m = $m should be ≥ 0."
    rtol ≥ 0 || @error "rtol = $rtol should be ≥ 0."
    atol ≥ 0 || @error "atol = $atol should be ≥ 0."
    maxit ≥ 0 || @error "maxit = $maxit should be ≥ 0."

    k = 0
    l₀ = norm(x)
    l₀ ≤ atol && return k, l₀, l₀  # lsol.m_active = [] falls to this as well

    τ = max(rtol*l₀, atol)

    n = length(x)
    m ≤ n || @error "m = $m should be > length(x) = $n."

    # Storing xold is needed only for the m ≠ 0 case, so executing the following block only
    # for m ≠ 0 would have saved 1 allocation.  However, putting the following block here
    # increase code readability, because doing update_lsol! right after norm_leq is an idiom.
    # m ≠ 0 is rarely used anyway.
    xold = Vector{T}(undef, n)
    xold .= x  # xold = x₀
    g!(x, xold)  # x is updated to x₁ = g(x₀)
    # Note we need at least two x's to perform the m ≠ 0 Anderson acceleration.

    if m == 0 # simple fixed-point iteration (no memory)
        while (k+=1) ≤ maxit
            l = norm(x, k, l₀)
            l ≤ τ && break

            # Evaluate g(xₖ).
            xold .= x  # xold = xₖ
            g!(x, xold)
        end
    else  # m ≠ 0
        # Pre-allocate all of the arrays we will need.  The goal is to allocate once and re-use
        # the storage during the iteration by operating in-place.
        f = Vector{T}(undef, n)
        ∆X = Matrix{T}(undef, n, m)  # columns: ∆x's
        ∆F = Matrix{T}(undef, n, m)  # columns: ∆f's
        Q = Matrix{T}(undef, n, m)  # space for QR factorization
        β = Vector{T}(undef, max(n,m))  # not m, in order to store RHS vector (f: length-n) and overwrite in-place via ldiv! (max length m)

        # Find x₁ = g(x₀): we need at least two x's to perform the Anderson acceleration.
        col = 1
        for i = 1:n
            f[i] = x[i] - xold[i]  # f(x₀) = g(x₀) - x₀
            ∆X[i,col] = f[i]  # ∆x₀ = x₁ - x₀
        end
        # Note that in the subsequent iteration, xₖ₊₁ is not necessarily g(xₖ), so
        # ∆xₖ = xₖ₊₁ - xₖ is not necessarily f(xₖ) = g(xₖ) - xₖ.

        # Perform subsequent iterations.
        #
        # (Assumption) The following quantities are known at the beginning of every iteration:
        # - xₖ
        # - f(xₖ₋₁)
        # - ∆xₖ₋₁
        # These quantities should be updated at the end of the for loop to satisfy the
        # assumption at every iteration.
        #
        # Using these quantities, in each iteration
        # - g(xₖ) is calculated,
        # - with which f(xₖ) = g(xₖ) - xₖ is calculated,
        # - with which ∆fₖ₋₁ = f(xₖ) - f(xₖ₋₁) is calculated.
        #
        # Then we know ∆xₖ₋₁ and ∆fₖ₋₁, which are needed for the Anderson acceleration.
        while (k+=1) ≤ maxit
            l = norm(x, k, l₀)
            l ≤ τ && break

            # Evaluate g(xₖ).
            xold .= x  # xold = xₖ
            g!(x, xold)

            # Prepare the least squares problem.
            for i = 1:n
                β[i] = x[i] - xold[i]  # temporary storage for f(xₖ) = g(xₖ) - xₖ
                ∆F[i,col] = β[i] - f[i]  # ∆fₖ₋₁ = f(xₖ) - f(xₖ₋₁)
                f[i] = β[i]  # f(xₖ)
            end

            # Solve the least squares problem.
            #
            # The code overwrites β with ∆F′ \ f.  (Note that β stores f at this point.)
            # Only first m′ entries of the solution β are meaningful, because they are the
            # entries of the least squares solution.
            #
            # In the following, QR′ and QR have different types, so we introduce separate
            # variables for type stability.
            if k < m
                # Construct subarrays to work in-place on a subset of the columns.
                Q′, ∆F′ = @views Q[:,1:k], ∆F[:,1:k]
                QR′ = qr!(copyto!(Q′, ∆F′), ColumnNorm())
                ldiv!(QR′, β)
            else
                QR = qr!(copyto!(Q, ∆F), ColumnNorm())
                ldiv!(QR, β)
            end

            # Replace the columns of ∆X (and ∆F in the next iteration) with the new data
            # in-place.  Rather than always appending the new data in the last column, cycle
            # through the m columns periodically.
            col = mod1(col+1, m)  # next column of ∆X and ∆F to update
            for i = 1:n
                # Update x with xₖ₊₁ = xₖ + fₖ - (∆X′ + ∆F′)*β′.
                # Note that x is already g(xₖ) = xₖ + f(xₖ)
                for j = 1:min(m,k)
                    x[i] -= (∆X[i,j] + ∆F[i,j]) * β[j]
                end

                ∆X[i,col] = x[i] - xold[i]  # ∆xₖ = xₖ₊₁ - xₖ
            end
        end  # for k = 1:maxit-1
    end  # if m == 0

    if k == maxit  # iteration terminated by consuming maxit steps
        l = norm(x, k, l₀)
        @warn "Anderson reached maxit = $maxit and didn't converge."
    end

    return k, l, l₀
end


text:
@discord_object struct Message
    id::Snowflake
    channel_id::Snowflake
    guild_id::Snowflake
    author::User
    member::GuildMember
    content::String
    timestamp::Union{String, DateTime}
    edited_timestamp::Union{String, DateTime}
    tts::Bool
    mention_everyone::Bool
    mentions::Vector{User}
    mention_roles::Vector{Snowflake}
    mention_channels::Vector{ChannelMention}
    attachments::Vector{Attachment}
    embeds::Vector{Embed}
    reactions::Vector{Reaction}
    nonce::Union{Int, String}
    pinned::Bool
    webhook_id::Snowflake
    type::MessageType.MessageTypeEnum
    activity::MessageActivity
    application::MessageApplication
    message_reference::MessageReference
    message_flags::Int
end


text:
module GainPatterns

import StatsBase.sample
import Base: +,-

export GainPattern, validgain, rotate!
export angular_error, angular_error_rel
export bearing_ls, bearing_cc, bearing_mle, bearing_sls, bearing_max
export bearing_ccn, bearing_half
export sample, rotate!, csv, normalize, normalize!, addsamples!, sampleGains
export samples_stds

_nullobs = 2147483647.

# angles is the list of angles
# meangains is the list of gains
# samples is a vector of gain vectors.
#  The mean of each of these makes up mean gains
type GainPattern

	angles::Vector{Float64}
	meangains::Vector{Float64}
	samples::Vector{Vector{Float64}}

	# Constructor where user passes in angles and mean gains
	# The samples array uses the single gain at each angle as that gain
	# TODO: check that these gains are valid
	# TODO: just don't insert that entry if invalid
	# WHY DO I MAKE IT A NAN??
	# Maybe it isn't totally useless info?
	function GainPattern(angles::Vector{Float64}, gains::Vector{Float64})

		# Error if the angles and gains vectors are of different lengths
		if length(angles) != length(gains)
			error("Number of angles must match number of gains")
		end

		# Create samples vectors from gains
		num_gains = length(gains)
		#samples = Array(Vector{Float64}, num_gains)
		samples = Array(Vector{Float64}, 0)
		new_angles = Array(Float64, 0)
		new_gains = Array(Float64, 0)
		for i = 1:num_gains
			if validgain(gains[i])
				push!(samples, [gains[i]])
				push!(new_angles, angles[i])
				push!(new_gains, gains[i])
				#samples[i] = [gains[i]]
			else
				#samples[i] = []
				#gains[i] = NaN
			end
		end
		return new(new_angles, new_gains, samples)
	end

	# Constructor taking in only angles and mean gains
	# Allows vectors of type real, and converts to float arrays
	function GainPattern{T1<:Real,T2<:Real}(angles::Vector{T1}, gains::Vector{T2})
		GainPattern(float(angles), float(gains))
	end

	# Allow user to input matrix
	# Makes no assumption about the validity of all the inputs
	function GainPattern{T1<:Real,T2<:Real}(angles::Vector{T1}, gains::Matrix{T2})

		# Determine number of angles and samples we are looking at
		num_angles, num_samples = size(gains)

		# Turn gains matrix into proper samples vector, then refine and return
		samples = Array(Vector{Float64}, num_angles)
		for i = 1:num_angles
			samples[i] =  vec(gains[i,:])
		end
		meangains, refined_samples = refine(samples)
		return new(angles, meangains, refined_samples)
	end

	# Create gain pattern from samples
	# First, we refine the samples (remove invalid entries, compute mean)
	# Then, just call the default constructor
	function GainPattern{T1<:Real,T2<:Real}(angles::Vector{T1}, samples::Vector{Vector{T2}})
		meangains, refined_samples = refine(samples)
		return new(angles, meangains, refined_samples)
	end

	# Create gain pattern from csv file
	# Assumes first column is angles
	# This cuts out any null observations..
	function GainPattern(datafile::AbstractString)

		# Read the data, creating the angles vector and an array of data
		alldata = readcsv(datafile)
		angles = float(alldata[:, 1])
		gains = alldata[:, 2:end]

		if size(gains, 2) == 1
			return GainPattern(angles, vec(gains))
		else
			return GainPattern(angles, gains)
		end

	end

end # type GainPattern

# Determes if gain value is valid
# In my research, I receive a value of 2147483647 if the gain is invalid
# If this is the value of the gain, I know it is not valid
# The user can override this to meet their own data needs
function validgain(gain::Real)
	#nullobs = 2147483647.
	gain == _nullobs ? false : true
end

# Eliminates any invalid gain values, and computes the mean
# Returns vector with means, and new vector with samples
function refine(samples::Vector{Vector{Float64}})

	# Calculate number of angles and create required data structures
	num_angles = length(samples)
	refined_samples = Array(Vector{Float64}, num_angles)
	meangains = Array(Float64, num_angles)

	# Loop through all sample vectors,
	#  eliminating invalid gain values and calculating the mean
	for i = 1:num_angles
		refined_samples[i] = Array(Float64, 0)
		rowavg = 0.0
		numvalid = 0
		num_samples = length(samples[i])
		for j = 1:num_samples
			tempgain = samples[i][j]
			push!(refined_samples[i], tempgain)
			if validgain(tempgain)
				rowavg += tempgain
				numvalid += 1
			end
		end
		meangains[i] = rowavg / numvalid
	end

	return meangains, refined_samples
end

# Subtracts gp2 from gp1
# Only works if they both have the same number of angles
# Assumes that the angles are the same
# TODO: Should I check that the angles are the same?
# If the samples are all of equal length, these are subtracted
# Otherwise, only the means are subtracted,
#  and a new samples array is initialized (with only the mean as a sample)
function -(gp1::GainPattern, gp2::GainPattern)
	num_angles1 = length(gp1.angles)
	num_angles2 = length(gp2.angles)

	if num_angles1 != num_angles2
		error("Number of angles must be the same!")
	end

	# Loop through samples, ensuring they are the same length
	same_sample_length = true
	num_angles = num_angles1
	gains = Array(Float64, num_angles)
	for i = 1:num_angles
		if length(gp1.samples[i]) != length(gp2.samples[i])
			same_sample_length = false
		end
		gains[i] = gp1.meangains[i] + gp2.meangains[i]
	end

	# If the samples are all the same length, add them
	# Otherwise, just create a new gain pattern with the calculated means
	samples = Array(Vector{Float64}, num_angles)
	if same_sample_length
		for i = 1:num_angles
			#samples[i] = gp1.samples[i] + gp2.samples[i]
			samples[i] = Array(Float64, 0)
			for j = 1:length(gp1.samples[1])
				if !validgain(gp1.samples[i][j]) || !validgain(gp2.samples[i][j])
					push!(samples[i], _nullobs)
				else
					push!(samples[i], gp1.samples[i][j] - gp2.samples[i][j])
				end
			end
		end
		gp_new = GainPattern(gp1.angles, samples)
	else
		gp_new = GainPattern(gp1.angles, gains)
	end

	return gp_new
end

# Addition function for gain patterns.
# Identical to subtraction, except for + sign in a few spots instead of -
# Leave invalid observations alone
function +(gp1::GainPattern, gp2::GainPattern)
	num_angles1 = length(gp1.angles)
	num_angles2 = length(gp2.angles)

	if num_angles1 != num_angles2
		error("Number of angles must be the same!")
	end

	# Loop through samples, ensuring they are the same length
	same_sample_length = true
	num_angles = num_angles1
	gains = Array(Float64, num_angles)
	for i = 1:num_angles
		if length(gp1.samples[i]) != length(gp2.samples[i])
			same_sample_length = false
		end
		gains[i] = gp1.meangains[i] + gp2.meangains[i]
	end

	# If the samples are all the same length, add them
	# Otherwise, just create a new gain pattern with the calculated means
	samples = Array(Vector{Float64}, num_angles)
	if same_sample_length
		for i = 1:num_angles
			#samples[i] = gp1.samples[i] + gp2.samples[i]
			samples[i] = Array(Float64, 0)
			for j = 1:length(gp1.samples[1])
				if !validgain(gp1.samples[i][j]) || !validgain(gp2.samples[i][j])
					push!(samples[i], _nullobs)
				else
					push!(samples[i], gp1.samples[i][j] + gp2.samples[i][j])
				end
			end
		end
		gp_new = GainPattern(gp1.angles, samples)
	else
		gp_new = GainPattern(gp1.angles, gains)
	end

	return gp_new
end

# Adding a constant term to a gain pattern
# Leave any invalid observations alone
function +(gp::GainPattern, c::Real)
	gp_new = deepcopy(gp)
	num_angles = length(gp.angles)
	for i = 1:num_angles
		gp_new.meangains[i] += c
		for j = 1:length(gp_new.samples[i])
			if validgain(gp_new.samples[i][j])
				gp_new.samples[i][j] += c
			end
		end
	end

	return gp_new
end

+(c::Real, gp::GainPattern) = +(gp, c)
-(gp::GainPattern, c::Real) = +(gp, -c)

# Appends the samples of gp2 to the samples vector of gp1
# Recalculates the mean gains of gp1
function addsamples!(gp1::GainPattern, gp2::GainPattern)
	@assert gp1.angles == gp2.angles
	num_angles = length(gp1.angles)
	for i = 1:num_angles
		for j = 1:length(gp2.samples[i])
			push!(gp1.samples[i], gp2.samples[i][j])
		end
	end
	meangains, refined_samples = refine(gp1.samples)
	gp1.meangains = meangains
end

# Returns a normalized distribution
# This normalization was shown in Graefenstein 2009
# Creates a copy of the gains then calls the in-place version
function normalize{T<:Real}(gains::Vector{T})
	new_gains = float(gains)	# create copy of gains, while ensuring float
	normalize!(new_gains)
	return new_gains
end

# Normalizes a vector by subtracting the mean and dividing by std dev
# Pretends invalid gains aren't even there
function normalize!(gains::Vector{Float64})

	temp_gains = copy(gains)
	len = length(temp_gains)
	for j = len:-1:1
		if !validgain(temp_gains[j])
			deleteat!(temp_gains, j)
		end
	end

	m = mean(temp_gains)
	s = std(temp_gains)

	for i = 1:len
		if validgain(gains[i])
			gains[i] = (gains[i] - m) / s
		end
	end
end

# Calculates mean and standard deviation of meangains
# Applies the normalization (x-m) / s
# Also applies normalization to each of the samples, but only valid ones
function normalize!(gp::GainPattern)
	m = mean(gp.meangains)
	s = std(gp.meangains)
	num_angles = length(gp.angles)

	for i = 1:num_angles
		gp.meangains[i] = (gp.meangains[i]-m)/s
		for j = 1:length(gp.samples[i])
			if validgain(gp.samples[i][j])
				gp.samples[i][j] = (gp.samples[i][j] - m) / s
			end
		end
	end
end

# Rotates the gain pattern by a specified number of degrees
# We must also call mod so the degrees stay in the proper range
function rotate!(gp::GainPattern, degrees::Real)
	for i = 1:length(gp.angles)
		gp.angles[i] = mod(gp.angles[i] + degrees, 360)
	end
end

###########################################################################
# ANGULAR ERROR
# Error returned is always a mangitude (positive quantity)
# Assumes both errors are between zero and 360
###########################################################################
function angular_error(angle1::Real, angle2::Real)
	angular_error(float(angle1), float(angle2))
end
function angular_error(angle1::Float64, angle2::Float64)
	err = angle2 - angle1
	if err > 180
		err = 360 - err
	elseif err < -180
		err += 360
	end
	if err < 0
		err *= -1
	end
	return err
end
export tangular_error
function tangular_error(angle1::Float64, angle2::Float64)
	temp = angle2 - angle1
	min( (abs(temp),angle2-angle1), (abs(360.0+temp),360+temp), (abs(temp-360),temp-360) )[2]
end
function angular_error{T<:Real}(angles1::Vector{T}, angles2::Vector{T})
	angular_error(float(angles1), float(angles2))
end
function angular_error(angles1::Vector{Float64}, angles2::Vector{Float64})
	@assert length(angles1) == length(angles2)
	err_len = length(angles1)
	err_arr = zeros(err_len)
	for i = 1:err_len
		err_arr[i] = angular_error(angles1[i], angles2[i])
	end
	return err_arr
end
function angular_error{T<:Real}(angles1::Vector{T}, angle2::Real)
	err_len = length(angles1)
	angles2 = angle2 * ones(err_len)
	return angular_error(float(angles1), angles2)
end
function angular_error{T<:Real}(angle1::Real, angles2::Vector{T})
	return angular_error(angles2, angle1)
end

# Finds the shortest distance between points, but preserves a sign
# say v2[i] = 30, v1[i] = 15 => ans[i] = 15
# say v2[i] = 15, v1[i] = 30 => ans[i] = -15
function angular_error_rel{T1<:Real, T2<:Real}(v1::Vector{T1}, v2::Vector{T2})
	@assert length(v1) == length(v2)
	err_len = length(v1)
	err_arr = zeros(err_len)
	for i = 1:err_len
		err_arr[i] = angular_error_rel(v1[i], v2[i])
	end
	return err_arr
end
function angular_error_rel(angle1::Float64, angle2::Float64)
	return tangular_error(angle1, angle2)
	#eps_val = 1e-6
	#ref = angular_error(angle1, angle2)
	#if abs(angle2 - angle1 - ref) < eps_val
	#	# we know it is positive
	#	ref = ref
	#elseif abs(360. + angle2 - angle1 - ref) < eps_val
	#	ref = ref
	#else
	#	ref = -ref
	#end
end
function angular_error_rel(angle1::Real, angle2::Real)
	angular_error_rel(float(angle1), float(angle2))
end

###########################################################################
# SAMPLING
###########################################################################
# Samples at a specified angle
# TODO: allow for angles not specified in pattern. Requires interpolation
function sample(gp::GainPattern, angle::Real)
	i = findfirst(gp.angles, angle)
	i == 0 ? error("Currently, you must sample from angle for which the gain pattern has samples") : nothing
	sample(gp.samples[i])
end

# Samples at a vector of specified angles
# These angles must all exist in the gain pattern's angle vector
# TODO: maybe allow for angles not specified in pattern?
function sample{T<:Real}(gp::GainPattern, angles::Vector{T})
	num_angles = length(angles)
	gains = Array(Float64, num_angles)
	for i = 1:num_angles
		idx = findfirst(gp.angles, angles[i])
		idx == 0 ? error("Currently, you must sample from angle for which the gain pattern has samples") : nothing
		gains[i] = sample(gp.samples[idx])
	end
	return gains
end

# Sample a gain for every angle in gp.
function sample(gp::GainPattern)
	sample(gp, gp.angles)
end

# Sample from a given distribution of gains
# Returns the gains sampled, and the angles at which these samples occurred
# TODO: Look into this, see if it is still working
function sampleGains(gains, step::Int64, bearing::Int64)
	noise = 0.25*randn(360)

	# antenna_angles are angles antenna is oriented at
	# sampled_angles is angle antenna observes in its measurements
	antenna_angles = [0:step:359]
	sampled_angles = mod(antenna_angles - bearing, 360)
	sampled_gains = gains[sampled_angles + 1] + noise[sampled_angles + 1]
	normalize!(sampled_gains)
	return sampled_gains, antenna_angles
end

###########################################################################
# INTERPOLATION
# Assumes ref_angles always given in 0-360 order
# First ref_angle will always be zero, but last will be less than 360
###########################################################################
# Finds the value (from ref_gains) at angle + shift
# Looks like linear interpolation to me
function interp_linear(angle, shift, ref_angles, ref_gains)
	# Add the shift to the desired angle, be sure to mod it
	angle = mod(angle+shift, 360)

	# Find the two ref_angles this is between
	# final i value will be between i, i+1
	# this is index of lower close value
	i = 1
	ref_len = length(ref_angles)
	while (i < ref_len) && (angle > ref_angles[i+1])
		i += 1
	end

	if i == ref_len
		# Between last and first ref_angles
		val = (360.0 - angle) * ref_gains[i]
		val += (angle - ref_angles[i]) * ref_gains[1]
		val = val / (360.0 - ref_angles[i])
	else
		val = (ref_angles[i+1] - angle) * ref_gains[i]
		val += (angle - ref_angles[i]) * ref_gains[i+1]
		val = val / (ref_angles[i+1] - ref_angles[i])
	end
	return val
end

function interp_nearest(angle, shift, ref_angles, ref_gains)

	# Determine angle and ref_angle indices it falls between
	angle = mod(angle+shift, 360)
	i = 1
	ref_len = length(ref_angles)
	while (i < ref_len) && (angle > ref_angles[i+1])
		i += 1
	end

	if i == ref_len
		# between last ref_angle and first ref_angle
		d1 = abs(ref_angles[i] - angle)
		d2 = abs(360 - angle)
		if d2 > d1
			val = ref_gains[i]
		else
			val = ref_gains[1]
		end
	else
		# Determine which you are closer too...
		d1 = abs(ref_angles[i] - angle)
		d2 = abs(ref_angles[i+1] - angle)
		if d2 > d2
			val = ref_gains[i]
		else
			val = ref_gains[i+1]
		end
	end
	return val
end

# Writes a gain pattern to a file
# TODO: check that the file exists
# TODO: Check that they add csv file extension if not, add it
# Currently only prints out meangains, but could theoretically do samples too
function csv(gp::GainPattern, filename::AbstractString="temp.csv"; samples=true)
	num_angles = length(gp.angles)
	outfile = open(filename, "w")

	if samples
		for i = 1:num_angles
			num_samples = length(gp.samples[i])
			write(outfile, "$(gp.angles[i])")
			for j = 1:num_samples
				write(outfile, ",$(gp.samples[i][j])")
			end
			write(outfile, "\n")
		end
	else
		for i = 1:num_angles
			write(outfile, "$(gp.angles[i]),$(gp.meangains[i])\n")
		end
	end

	close(outfile)
end

# Get standard deviation of the each sample
# Assumes the meangains value is corrrect (and exists)
function samples_stds(gp::GainPattern)
	num_angles = length(gp.angles)
	stds = zeros(num_angles)
	for i = 1:num_angles
		num_samples = length(gp.samples[i])
		sample_var = 0.0
		angle_mean = gp.meangains[i]
		n = 0.0
		for j = 1:num_samples
			temp = gp.samples[i][j]
			if validgain(temp)
				sample_var += (temp - angle_mean)^2
				n += 1.0
			end
		end
		stds[i] = sample_var / n
	end
	return sqrt(stds)
end

###########################################################################
# EXTREMA
# I don't remember what I used this for, but it appears unused
###########################################################################
# Finds minimum and maximum gains in gain pattern's samples, in one sweep
# Checks if each gain is valid
# I don't know how I feel about that...
function Base.extrema(gp::GainPattern)
	num_angles = length(gp.angles)
	min_sample = Inf
	max_sample = -Inf
	for i = 1:num_angles
		num_samples = length(gp.samples[i])
		for j = 1:num_samples
			temp = gp.samples[i][j]
			if validgain(temp)
				max_sample = max(temp, max_sample)
				min_sample = min(temp, min_sample)
			end
		end
	end
	return min_sample, max_sample
end

# Returns the maximum or minimum sample in the gain pattern's samples
Base.maximum(gp::GainPattern) = extrema(gp)[2]
Base.minimum(gp::GainPattern) = extrema(gp)[1]

# Methods to compute bearing
include("bearing/cross_correlation.jl")
include("bearing/half.jl")
include("bearing/mle.jl")
include("bearing/least_squares.jl")
include("bearing/scaled_least_squares.jl")
include("bearing/max.jl")

end # module


text:
using CompScienceMeshes, BEAST
Γ = readmesh(joinpath(@__DIR__,"sphere2.in"))

X = raviartthomas(Γ)
Y = buffachristiansen(Γ)

Δt, Nt = 0.6, 200
T = timebasisshiftedlagrange(Δt, Nt, 2)
δ = timebasisdelta(Δt, Nt)

V = X ⊗ T
W = Y ⊗ δ

duration = 20 * Δt
delay = 1.5 * duration
amplitude = 1.0
gaussian = creategaussian(duration, delay, amplitude)

direction, polarisation = ẑ, x̂
E = BEAST.planewave(polarisation, direction, gaussian, 1.0)
H = direction × E

K = TDMaxwell3D.doublelayer(speedoflight=1.0)
I = Identity()
N = NCross()

@hilbertspace k
@hilbertspace j
mfie = @discretise (0.5(N⊗I) + 1.0K)[k,j] == -1.0H[k] k∈W j∈V
xmfie = solve(mfie)

Xmfie, Δω, ω0 = fouriertransform(xmfie, Δt, 0.0, 2)
ω = collect(ω0 .+ (0:Nt-1)*Δω)
_, i1 = findmin(abs.(ω .- 1.0))

ω1 = ω[i1]
um = Xmfie[:,i1] / fouriertransform(gaussian)(ω1)


text:
# g = grammar_from_string(
#     """
#     S --> NP VP
#     NP --> D N | dog
#     VP --> V NP | V PP
#     PP --> P N
#     D --> the | prep
#     N --> dog
#     V --> paints
#     P --> prep
#     """
# );
g = grammar_from_string( # N --> T1 and so on don't work ; T1 treated as a terminal!
    """
    N --> hand | size | ADJ T6 | V T9
    ADJ --> great | N T7 | V T8
    V --> drench | ADJ T10
    ADV --> very | ADJ T11
    T1 --> hand
    T2 --> size
    T3 --> great
    T4 --> drench
    T5 --> very
    T6 --> +ness
    T7 --> +y
    T8 --> +ed
    T9 --> +er
    T10 --> +ify
    T11 --> +ly
    """
);
g.start_categories = g.categories[1:4]

test_str = "hand +y"

function parseString(grammar, str) #get all parses for string
    scores = parse(grammar, split(str))
    g2 = add_score(grammar, :enum_forest, enum_forest_score)
    score = parse(g2, split(str))["S"]
    trees = tree_structs(g, score.enum_forest)
    return trees
end

# test_str = "the dog paints prep dog"

function sampleTree(grammar, str) #sample a parse from possble parses of string
    g3 = add_forest_score(add_random_prob_score(grammar, :count), :prob)
    forest = parse(g3, split(str))["S"].forest;
    rand(forest) |>            # get random tree iterator
    collect |>                 # collect rule applications of a random tree
    t -> tree_struct(g3, t)
end

#sampleTree(g, "the dog paints prep dog")

#function forwardSample(grammar, dist)

# v = [5, 10, 15]
# r = ChineseRest(v)
#println(r)
#println(g)


text:
module Galerkin

export
    galerkin_projection,
    galerkin_projection!,
    galerkin_matrix,
    galerkin_matrix!,
    galerkin_tensor,
    galerkin_tensor!

using ..BandedTensors
using ..BSplines
using ..DifferentialOps
using ..Recombinations: num_constraints
using ..Collocation: collocation_matrix  # for Documenter

using BandedMatrices
using FastGaussQuadrature: gausslegendre
using LinearAlgebra: Hermitian, ⋅
using SparseArrays
using StaticArrays

const DerivativeCombination{N} = Tuple{Vararg{Derivative,N}}

include("quadratures.jl")
include("projection.jl")  # galerkin_projection
include("linear.jl")      # galerkin_matrix
include("quadratic.jl")   # galerkin_tensor

end

max_stars_repo_name: jipolanco/BasisSplines

max_issues_repo_name: jipolanco/BasisSplines

max_forks_repo_name: jipolanco/BasisSplines


text:
"""
   probability_check(pomdp::POMDP)

Checks if the transition and observation function of the discrete `pomdp`
have probability mass that sums up to unity for all state-action pairs.
"""
function probability_check(pomdp::POMDP)
    obs_prob_consistency_check(pomdp)
    trans_prob_consistency_check(pomdp)
end

"""
    obs_prob_consistency_check(pomdp::POMDP)

Checks if the observation function of the discrete `pomdp`
has probability mass that sums up to unity for all state-action pairs.
"""
function obs_prob_consistency_check(pomdp::POMDP)
    # initalize space
    sspace = states(pomdp)
    aspace = actions(pomdp)
    ospace = observations(pomdp)
    # iterate through all s-a pairs
    for s in iterator(sspace)
        for a in iterator(aspace)
            obs = observation(pomdp, a, s)
            p = 0.0
            for o in iterator(ospace)
                p += pdf(obs, o)
            end
            @assert isapprox(p, 1.0) "Observation probability does not sum to unity for state: ", s, " action: ", a
        end
    end
end

"""
    trans_prob_consistency_check(pomdp::Union{MDP, POMDP})

Checks if the transition function of the discrete problem
has probability mass that sums up to unity for all state-action pairs.
"""
function trans_prob_consistency_check(pomdp::Union{MDP, POMDP})
    # initalize space
    sspace = states(pomdp)
    aspace = actions(pomdp)
    # iterate through all s-a pairs
    for s in iterator(sspace)
        for a in iterator(aspace)
            tran = transition(pomdp, s, a)
            p = 0.0
            for sp in iterator(sspace)
                p += pdf(tran, sp)
            end
            @assert isapprox(p, 1.0) "Transition probability does not sum to unity for state: ", s, " action: ", a
        end
    end
end


text:
#=
Given the mapping a = 1, b = 2, ... z = 26, and an encoded message, count the number of ways it can be decoded.

For example, the message '111' would give 3, since it could be decoded as 'aaa', 'ka', and 'ak'.

You can assume that the messages are decodable. For example, '001' is not allowed.
=#

function iscode(code)
    if typeof(code) == String
        code = parse(Int, code)
    end
    if code >= 1 && code <= 26
        return true
    else
        return false
    end
end

function count_ways(code)
    if typeof(code) == String
        code = parse(Int, code)
    end
    length_code = length(string(code))
    count = 0
    if length_code == 1
        count = 1
    elseif length_code == 2 && iscode(code)
        count = 2
    elseif length_code == 2 && iscode(code) == false
        count = 1
    else
        if iscode(string(code)[1:2])
            count = count_ways(string(code)[3:end]) + count
        end
        count = count_ways(string(code)[2:end]) + count
    end
    return count
end

max_stars_repo_name: DominiqueCaron/daily-coding-problem

max_issues_repo_name: DominiqueCaron/daily-coding-problem

max_forks_repo_name: DominiqueCaron/daily-coding-problem


text:
"""
    project(u::AbstractVector, v::AbstractVector) -> AbstractVector

Project vector `u` onto vector `v`.

# Examples

```jldoctest
julia> project([1, 1], [1, 0])
2-element Base.Vector{Float64}:
 1.0
 0.0

julia> project([5, 5], [1, 0])
2-element Base.Vector{Float64}:
 5.0
 0.0

julia> project([5, -5], [0, 1])
2-element Base.Vector{Float64}:
 -0.0
 -5.0
```

"""
function project(u::AbstractVector, v::AbstractVector)
    return (u ⋅ v) / (v ⋅ v) * v
end

"""
    project(point::AbstractVector, line::AbstractLine) -> StaticArrays.SVector

Project a point onto a line.

# Examples

```jldoctest
julia> project([1, 1], Line([0, 0], [1, 0]))
2-element StaticArrays.SVector{2, Float64} with indices SOneTo(2):
 1.0
 0.0

julia> project([5, -1], Line([0, 0], [1, 0]))
2-element StaticArrays.SVector{2, Float64} with indices SOneTo(2):
 5.0
 0.0

julia> project([1, 0], Line([0, 0], [1, 1]))
2-element StaticArrays.SVector{2, Float64} with indices SOneTo(2):
 0.5
 0.5

julia> point = project([1, 0, 0], Line([0, 0, 0], [1, 1, 1]));

julia> round.(point, digits=3)
3-element StaticArrays.SVector{3, Float64} with indices SOneTo(3):
 0.333
 0.333
 0.333
```

"""
function project(point::AbstractVector, line::AbstractLine)

    # Vector from the point on the line to the point in space.
    vector_to_point = Vector(line.point, point)

    # Vector projected onto the line.
    vector_projected = project(vector_to_point, line.direction)

    return line.point + vector_projected
end

"""
    project(point::AbstractVector, plane::AbstractPlane) -> StaticArrays.SVector

Project a point onto a plane.

# Examples

```jldoctest
julia> project([0, 0, 5], Plane([0, 0, 0], [0, 0, 1]))
3-element StaticArrays.SVector{3, Float64} with indices SOneTo(3):
 0.0
 0.0
 0.0

julia> plane = Plane([1,2,3], [1, 3, -2]);

julia> point_projected = project([5, 1, 3], plane);

julia> round.(point_projected, digits=3)
3-element StaticArrays.SVector{3, Float64} with indices SOneTo(3):
 4.929
 0.786
 3.143
```

"""
function project(point::AbstractVector, plane::AbstractPlane)

    # Vector from the point in space to the point on the plane.
    vector_to_plane_point = Vector(point, plane.point)

    # Perpendicular vector from the point in space to the plane.
    vector_to_plane = project(vector_to_plane_point, plane.normal)

    return point + vector_to_plane
end


text:
using SIAN

println("Setting up the problem")

ode = @ODEmodel(
            x1'(t) = -beta * x1(t) * x4(t) - d * x1(t) + s,
            x2'(t) = beta * q1 * x1(t) * x4(t) - k1 * x2(t) - mu1 * x2(t),
            x3'(t) = beta * q2 * x1(t) * x4(t) + k1 * x2(t) - mu2 * x3(t),
            x4'(t) = -c * x4(t) + k2 * x3(t),
            y1(t) = x1(t),
            y2(t) = x4(t)
)

res = identifiability_ode(ode, get_parameters(ode); p=0.99, p_mod=0, nthrds=1)

println(res)

max_stars_repo_name: iliailmer/SIAN-Julia

max_issues_repo_name: iliailmer/SIAN-Julia

max_forks_repo_name: iliailmer/SIAN-Julia


text:
# This file includes code that was formerly a part of Julia. License is MIT: http://julialang.org/license

## required core functionality ##

lastindex(s::ASCIIString) = length(s.data)
getindex(s::ASCIIString, i::Int) = (x=s.data[i]; ifelse(x < 0x80, Char(x), '\ufffd'))

codeunit(s::ASCIIString) = UInt8
ncodeunits(s::ASCIIString) = length(s.data)

if isdefined(Base, :iterate)
    import Base: iterate
    function iterate(s::ASCIIString, i::Int = firstindex(s))
        i > ncodeunits(s) && return nothing
        return next(s, i)
    end
end

## overload methods for efficiency ##

bytestring(s::ASCIIString) = s

sizeof(s::ASCIIString) = sizeof(s.data)

getindex(s::ASCIIString, r::Vector) = ASCIIString(getindex(s.data,r))
getindex(s::ASCIIString, r::UnitRange{Int}) = ASCIIString(getindex(s.data,r))
getindex(s::ASCIIString, indx::AbstractVector{Int}) = ASCIIString(s.data[indx])
function search(s::ASCIIString, c::Char, i::Integer)
    i == sizeof(s) + 1 && return 0
    (i < 1 || i > sizeof(s)) && throw(BoundsError(s, i))
    return c < Char(0x80) ? search(s.data,c%UInt8,i) : 0
end
rsearch(s::ASCIIString, c::Char, i::Integer) = c < Char(0x80) ? rsearch(s.data,c%UInt8,i) : 0

function string(c::ASCIIString...)
    if length(c) == 1
        return c[1]
    end
    n = 0
    for s in c
        n += length(s.data)
    end
    v = Vector{UInt8}(undef, n)
    o = 1
    for s in c
        ls = length(s.data)
        unsafe_copy!(v, o, s.data, 1, ls)
        o += ls
    end
    ASCIIString(v)
end

function ucfirst(s::ASCIIString)
    if !isempty(s) && 'a' <= s[1] <= 'z'
        t = ASCIIString(copy(s.data))
        t.data[1] -= 32
        return t
    end
    return s
end
function lcfirst(s::ASCIIString)
    if !isempty(s) && 'A' <= s[1] <= 'Z'
        t = ASCIIString(copy(s.data))
        t.data[1] += 32
        return t
    end
    return s
end

function uppercase(s::ASCIIString)
    d = s.data
    for i = 1:length(d)
        if 'a' <= Char(d[i]) <= 'z'
            td = copy(d)
            for j = i:length(td)
                if 'a' <= Char(td[j]) <= 'z'
                    td[j] -= 32
                end
            end
            return ASCIIString(td)
        end
    end
    return s
end
function lowercase(s::ASCIIString)
    d = s.data
    for i = 1:length(d)
        if 'A' <= Char(d[i]) <= 'Z'
            td = copy(d)
            for j = i:length(td)
                if 'A' <= Char(td[j]) <= 'Z'
                    td[j] += 32
                end
            end
            return ASCIIString(td)
        end
    end
    return s
end

reverse(s::ASCIIString) = ASCIIString(reverse(s.data))

## outputing ASCII strings ##

write(io::IO, s::ASCIIString) = write(io, s.data)

## transcoding to ASCII ##

ascii(x) = convert(ASCIIString, x)
convert(::Type{ASCIIString}, s::ASCIIString) = s
convert(::Type{ASCIIString}, s::String) = ascii(codeunits(s))
convert(::Type{ASCIIString}, s::UTF8String) = ascii(s.data)
convert(::Type{ASCIIString}, a::Vector{UInt8}) = begin
    isvalid(ASCIIString,a) || throw(ArgumentError("invalid ASCII sequence"))
    return ASCIIString(a)
end
if isdefined(Base, :codeunits)
    convert(::Type{ASCIIString}, a::Base.CodeUnits{UInt8,String}) = convert(ASCIIString, Vector{UInt8}(a))
end

ascii(p::Ptr{UInt8}) =
    ascii(p, p == C_NULL ? Csize_t(0) : ccall(:strlen, Csize_t, (Ptr{UInt8},), p))
function ascii(p::Ptr{UInt8}, len::Integer)
    p == C_NULL && throw(ArgumentError("cannot convert NULL to string"))
    ary = ccall(:jl_pchar_to_array, Vector{UInt8},
                (Ptr{UInt8}, Csize_t), p, len)
    isvalid(ASCIIString, ary) || throw(ArgumentError("invalid ASCII sequence"))
    ASCIIString(ary)
end

function convert(::Type{ASCIIString}, a::Array{UInt8,1}, invalids_as::ASCIIString)
    l = length(a)
    idx = 1
    iscopy = false
    while idx <= l
        (a[idx] < 0x80) && (idx +=1; continue)
        !iscopy && (a = copy(a); iscopy = true)
        endn = idx
        while endn <= l
            (a[endn] < 0x80) && break
            endn += 1
        end
        (endn > idx) && (endn -= 1)
        splice!(a, idx:endn, invalids_as.data)
        l = length(a)
    end
    convert(ASCIIString, a)
end
convert(::Type{ASCIIString}, a::Array{UInt8,1}, invalids_as::AbstractString) =
    convert(ASCIIString, a, ascii(invalids_as))
convert(::Type{ASCIIString}, s::AbstractString) = ascii(bytestring(s))


text:
@testset "Operators" begin
    println("Testing operators...")

    @testset "2D operators" begin
        println("  Testing 2D operators...")
        Nx, Ny, Nz = 32, 16, 8
        Lx, Ly, Lz = 100, 100, 100

        arch = CPU()
        grid = RegularCartesianGrid((Nx, Ny, Nz), (Lx, Ly, Lz))

        Hx, Hy, Hz = grid.Hx, grid.Hy, grid.Hz
        Tx, Ty, Tz = grid.Tx, grid.Ty, grid.Tz

        A3 = OffsetArray(zeros(Tx, Ty, Tz), 1-Hx:Nx+Hx, 1-Hy:Ny+Hy, 1-Hz:Nz+Hz)
        @. @views A3[1:Nx, 1:Ny, 1:Nz] = rand()
        fill_halo_regions!(A3, HorizontallyPeriodicBCs(), arch, grid)

        # A yz-slice with Nx==1.
        A2yz = OffsetArray(zeros(1+2Hx, Ty, Tz), 1-Hx:1+Hx, 1-Hy:Ny+Hy, 1-Hz:Nz+Hz)
        grid_yz = RegularCartesianGrid((1, Ny, Nz), (Lx, Ly, Lz))

        # Manually fill in halos for the slice.
        A2yz[0:2, 0:Ny+1, 1:Nz] .= A3[1:1, 0:Ny+1, 1:Nz]
        A2yz[:, :, 0] .= A2yz[:, :, 1]
        A2yz[:, :, Nz+1] .= A2yz[:, :, Nz]

        # An xz-slice with Ny==1.
        A2xz = OffsetArray(zeros(Tx, 1+2Hy, Tz), 1-Hx:Nx+Hx, 1-Hy:1+Hy, 1-Hz:Nz+Hz)
        grid_xz = RegularCartesianGrid((Nx, 1, Nz), (Lx, Ly, Lz))

        # Manually fill in halos for the slice.
        A2xz[0:Nx+1, 0:2, 1:Nz] .= A3[0:Nx+1, 1:1, 1:Nz]
        A2xz[:, :, 0] .= A2xz[:, :, 1]
        A2xz[:, :, Nz+1] .= A2xz[:, :, Nz]

        test_indices_3d = [(4, 5, 5), (21, 11, 4), (16, 8, 4),  (30, 12, 3), (11, 3, 6), # Interior
                           (2, 10, 4), (31, 5, 6), (10, 2, 4), (17, 15, 5), (17, 10, 2), (23, 5, 7),  # Borderlands
                           (1, 5, 5), (32, 10, 3), (16, 1, 4), (16, 16, 4), (16, 8, 1), (16, 8, 8),  # Edges
                           (1, 1, 1), (32, 16, 8)] # Corners

        test_indices_2d_yz = [(1, 1, 1), (1, 1, 2), (1, 2, 1), (1, 2, 2),
                              (1, 1, 5), (1, 5, 1), (1, 5, 5), (1, 11, 4),
                              (1, 15, 7), (1, 15, 8), (1, 16, 7), (1, 16, 8)]

        test_indices_2d_xz = [(1, 1, 1), (1, 1, 2), (2, 1, 1), (2, 1, 2),
                              (1, 1, 5), (5, 1, 1), (5, 1, 5), (17, 1, 4),
                              (31, 1, 7), (31, 1, 8), (32, 1, 7), (32, 1, 8)]

        for idx in test_indices_2d_yz
            @test δx_f2c(grid_yz, A2yz, idx...) ≈ 0
            @test δx_c2f(grid_yz, A2yz, idx...) ≈ 0
            @test δy_f2c(grid_yz, A2yz, idx...) ≈ δy_f2c(grid_yz, A3, idx...)
            @test δy_c2f(grid_yz, A2yz, idx...) ≈ δy_c2f(grid_yz, A3, idx...)
            @test δz_f2c(grid_yz, A2yz, idx...) ≈ δz_f2c(grid_yz, A3, idx...)
            @test δz_c2f(grid_yz, A2yz, idx...) ≈ δz_c2f(grid_yz, A3, idx...)
        end

        for idx in test_indices_2d_xz
            @test δx_f2c(grid_xz, A2xz, idx...) ≈ δx_f2c(grid_xz, A3, idx...)
            @test δx_c2f(grid_xz, A2xz, idx...) ≈ δx_c2f(grid_xz, A3, idx...)
            @test δy_f2c(grid_xz, A2xz, idx...) ≈ 0
            @test δy_c2f(grid_xz, A2xz, idx...) ≈ 0
            @test δz_f2c(grid_xz, A2xz, idx...) ≈ δz_f2c(grid_xz, A3, idx...)
            @test δz_c2f(grid_xz, A2xz, idx...) ≈ δz_c2f(grid_xz, A3, idx...)
        end
    end
end


text:
### Handle differences between versions. --------------------------------------------------------

if VERSION < v"0.4.0-dev+2418"
    # returns the index of the previous element for which the function returns true, or zero if it never does
    function findprev(testf::Function, A, start)
        for i = start:-1:1
            testf(A[i]) && return i
        end
        0
    end
    findlast(testf::Function, A) = findprev(testf, A, length(A))
end

if VERSION < v"0.4.0-dev+4393"
    # Return a relative filepath to path either from the current directory or from an optional start directory.
    # This is a path computation: the filesystem is not accessed to confirm the existence or nature of path or startpath.
    function relpath(path::ByteString, startpath::ByteString = ".")
        isempty(path) && throw(ArgumentError("`path` must be specified"))
        isempty(startpath) && throw(ArgumentError("`startpath` must be specified"))
        curdir = "."
        pardir = ".."
        path == startpath && return curdir

        path_arr  = split(abspath(path),      Base.path_separator_re)
        start_arr = split(abspath(startpath), Base.path_separator_re)

        i = 0
        while i < min(length(path_arr), length(start_arr))
            i += 1
            if path_arr[i] != start_arr[i]
                i -= 1
                break
            end
        end

        pathpart = join(path_arr[i+1:findlast(x -> !isempty(x), path_arr)], Base.path_separator)
        prefix_num = findlast(x -> !isempty(x), start_arr) - i - 1
        if prefix_num >= 0
            prefix = pardir * Base.path_separator
            relpath_ = isempty(pathpart)     ?
                (prefix^prefix_num) * pardir :
                (prefix^prefix_num) * pardir * Base.path_separator * pathpart
        else
            relpath_ = pathpart
        end
        return isempty(relpath_) ? curdir :  relpath_
    end
end

if VERSION < v"0.4.0-dev+4499"
    function cptree(src::AbstractString, dst::AbstractString; remove_destination::Bool=false,
                                                                 follow_symlinks::Bool=false)
        isdir(src) || throw(ArgumentError("'$src' is not a directory. Use `cp(src, dst)`"))
        if ispath(dst)
            if remove_destination
                rm(dst; recursive=true)
            else
                throw(ArgumentError(string("'$dst' exists. `remove_destination=true` ",
                                           "is required to remove '$dst' before copying.")))
            end
        end
        mkdir(dst)
        for name in readdir(src)
            srcname = joinpath(src, name)
            if !follow_symlinks && islink(srcname)
                symlink(readlink(srcname), joinpath(dst, name))
            elseif isdir(srcname)
                cptree(srcname, joinpath(dst, name); remove_destination=remove_destination,
                                                     follow_symlinks=follow_symlinks)
            else
                Base.FS.sendfile(srcname, joinpath(dst, name))
            end
        end
    end

    function cp(src::AbstractString, dst::AbstractString; remove_destination::Bool=false,
                                                             follow_symlinks::Bool=false)
        if ispath(dst)
            if remove_destination
                rm(dst; recursive=true)
            else
                throw(ArgumentError(string("'$dst' exists. `remove_destination=true` ",
                                           "is required to remove '$dst' before copying.")))
            end
        end
        if !follow_symlinks && islink(src)
            symlink(readlink(src), dst)
        elseif isdir(src)
            cptree(src, dst; remove_destination=remove_destination, follow_symlinks=follow_symlinks)
        else
            Base.FS.sendfile(src, dst)
        end
    end
end

if endswith(functionloc(isgeneric)[1], "deprecated.jl")
    _isgeneric(x) = true
else
    _isgeneric(x) = isgeneric(x)
end

lsdfield(x :: Function, f)    = _isgeneric(x) ? lsdfield(methods(x), f) : lsdfield(x.code, f)
lsdfield(x :: Method, f)      = lsdfield(x.func, f)
lsdfield(x :: MethodTable, f) = lsdfield(x.defs, f)

lsdfield(x :: LambdaStaticData, f) = getfield(x, f)


text:
"""
    InvWarpedView(img, tinv, [indices]; kwargs...) -> wv
    InvWarpedView(inner_view, tinv) -> wv

Create a view of `img` that lazily transforms any given index `I`
passed to `wv[I]` so that `wv[I] == img[inv(tinv)(I)]`.

Except for the lazy evaluation, the following two lines are expected to be equivalent:

```julia
warp(img, inv(tform), [indices]; kwargs...)
invwarpedview(img, tform, [indices]; kwargs...)
```

The conceptual difference to [`WarpedView`](@ref) is that
`InvWarpedView` is intended to be used when reasoning about the
image is more convenient that reasoning about the indices.
Furthermore, `InvWarpedView` allows simple nesting of
transformations, in which case the transformations will be
composed into a single one.

For detailed explaination of warp, associated arguments and parameters,
please refer to [`warp`](@ref).
"""
struct InvWarpedView{T,N,A,F,I,FI<:Transformation,E} <: AbstractArray{T,N}
    inner::WarpedView{T,N,A,F,I,E}
    inverse::FI
end

function InvWarpedView(inner::WarpedView{T,N,TA,F,I,E}) where {T,N,TA,F,I,E}
    tinv = inv(inner.transform)
    InvWarpedView{T,N,TA,F,I,typeof(tinv),E}(inner, tinv)
end

function InvWarpedView(A::AbstractArray, tinv::Transformation, inds::Tuple = autorange(A, tinv); kwargs...)
    inner = WarpedView(A, inv(tinv), inds; kwargs...)
    InvWarpedView(inner, tinv)
end

function InvWarpedView(inner::InvWarpedView, outer_tinv::Transformation)
    tinv = compose(outer_tinv, inner.inverse)
    InvWarpedView(parent(inner), tinv)
end

function InvWarpedView(inner::InvWarpedView, outer_tinv::Transformation, inds::Tuple)
    tinv = compose(outer_tinv, inner.inverse)
    InvWarpedView(parent(inner), tinv, inds)
end

@inline Base.parent(A::InvWarpedView) = parent(A.inner)
@inline Base.axes(A::InvWarpedView) = axes(A.inner)
@inline Base.size(A::InvWarpedView) = size(A.inner)

Base.@propagate_inbounds Base.getindex(A::InvWarpedView{T,N}, I::Vararg{Int,N}) where {T,N} = A.inner[I...]

function Base.showarg(io::IO, A::InvWarpedView, toplevel)
    print(io, "InvWarpedView(")
    Base.showarg(io, parent(A), false)
    print(io, ", ")
    print(io, A.inverse)
    if toplevel
        print(io, ") with eltype ", eltype(parent(A)))
    else
        print(io, ')')
    end
end


text:
# for this type of union, the compiler will actually generate the necessary methods
# but it's also type stable either way
_MdTypes = Union{BlockDescription,LikelihoodRatioTest,OptSummary,VarCorr,MixedModel}
Base.show(mime::MIME, x::_MdTypes) = show(Base.stdout, mime, x)

Base.show(io::IO, ::MIME"text/markdown", x::_MdTypes) = show(io, Markdown.MD(_markdown(x)))
# let's not discuss why we need show above and println below,
# nor what happens if we try display instead :)
Base.show(io::IO, ::MIME"text/html", x::_MdTypes) = println(io, Markdown.html(_markdown(x)))
# print and println because Julia already adds a newline line
Base.show(io::IO, ::MIME"text/latex", x::_MdTypes) = print(io, Markdown.latex(_markdown(x)))
function Base.show(io::IO, ::MIME"text/xelatex", x::_MdTypes)
    return print(io, Markdown.latex(_markdown(x)))
end

# not sure why this escaping doesn't work automatically
# FIXME: find out a way to get the stdlib to do this
function Base.show(io::IO, ::MIME"text/html", x::OptSummary)
    out = Markdown.html(_markdown(x))
    out = replace(out, r"`([^[:space:]]*)`" => s"\1")
    out = replace(out, r"\*\*(.*?)\*\*" => s"\1")
    return println(io, out)
end

function Base.show(io::IO, ::MIME"text/latex", x::OptSummary)
    out = Markdown.latex(_markdown(x))
    out = replace(out, r"`([^[:space:]]*)`" => s"\\texttt{\1}")
    out = replace(out, r"\*\*(.*?)\*\*" => s"\\textbf{\1}")
    return print(io, out)
end

function Base.show(io::IO, ::MIME"text/latex", x::MixedModel)
    la = Markdown.latex(_markdown(x))
    # take advantage of subscripting
    # including preceding & prevents capturing coefficients
    la = replace(la, r"& σ\\_([[:alnum:]]*) " => s"& $\\sigma_\\text{\1}$ ")
    return print(io, la)
end

function Base.show(io::IO, ::MIME"text/latex", x::LikelihoodRatioTest)
    la = Markdown.latex(_markdown(x))
    # take advantage of subscripting
    # including preceding & prevents capturing coefficients
    la = replace(la, r"χ²" => s"$\\chi^2$")
    return print(io, la)
end

function _markdown(b::BlockDescription)
    ncols = length(b.blknms)
    align = repeat([:l], ncols + 1)
    newrow = ["rows"; [bn for bn in b.blknms]]
    rows = [newrow]

    for (i, r) in enumerate(b.blkrows)
        newrow = [string(r)]
        for j in 1:i
            push!(newrow, "$(b.ALtypes[i, j])")
        end
        i < ncols && append!(newrow, repeat([""], ncols - i))
        push!(rows, newrow)
    end

    tbl = Markdown.Table(rows, align)
    return tbl
end

function _markdown(lrt::LikelihoodRatioTest)
    Δdf = lrt.tests.dofdiff
    Δdev = lrt.tests.deviancediff

    nr = length(lrt.formulas)
    outrows = Vector{Vector{String}}(undef, nr + 1)

    outrows[1] = ["", "model-dof", "deviance", "χ²", "χ²-dof", "P(>χ²)"] # colnms

    outrows[2] = [
        string(lrt.formulas[1]),
        string(lrt.dof[1]),
        string(round(Int, lrt.deviance[1])),
        " ",
        " ",
        " ",
    ]

    for i in 2:nr
        outrows[i + 1] = [
            string(lrt.formulas[i]),
            string(lrt.dof[i]),
            string(round(Int, lrt.deviance[i])),
            string(round(Int, Δdev[i - 1])),
            string(Δdf[i - 1]),
            string(StatsBase.PValue(lrt.pvalues[i - 1])),
        ]
    end

    tbl = Markdown.Table(outrows, [:l, :r, :r, :r, :r, :l])
    return tbl
end

_dname(::GeneralizedLinearMixedModel) = "Dispersion"
_dname(::LinearMixedModel) = "Residual"

function _markdown(m::MixedModel)
    if m.optsum.feval < 0
        @warn("Model has not been fit: results will be nonsense")
    end
    n, p, q, k = size(m)
    REML = m.optsum.REML
    nrecols = length(fnames(m))

    digits = 4

    co = coef(m)
    se = stderror(m)
    z = co ./ se
    p = ccdf.(Chisq(1), abs2.(z))

    σvec = vcat(collect.(values.(values(m.σs)))...)
    σwidth = _printdigits(σvec)

    newrow = ["", "Est.", "SE", "z", "p"]
    align = [:l, :r, :r, :r, :r]

    for rr in fnames(m)
        push!(newrow, "σ_$(rr)")
        push!(align, :r)
    end

    rows = [newrow]

    for (i, bname) in enumerate(coefnames(m))
        newrow = [
            bname,
            Ryu.writefixed(co[i], digits),
            Ryu.writefixed(se[i], digits),
            sprint(show, StatsBase.TestStat(z[i])),
            sprint(show, StatsBase.PValue(p[i])),
        ]
        bname = Symbol(bname)

        for (j, sig) in enumerate(m.σs)
            if bname in keys(sig)
                push!(newrow, Ryu.writefixed(getproperty(sig, bname), digits))
            else
                push!(newrow, " ")
            end
        end
        push!(rows, newrow)
    end

    re_without_fe = setdiff(
        mapfoldl(x -> Set(getproperty(x, :cnames)), ∪, m.reterms), coefnames(m)
    )

    for bname in re_without_fe
        newrow = [bname, "", "", "", ""]
        bname = Symbol(bname)

        for (j, sig) in enumerate(m.σs)
            if bname in keys(sig)
                push!(newrow, Ryu.writefixed(getproperty(sig, bname), digits))
            else
                push!(newrow, " ")
            end
        end
        push!(rows, newrow)
    end

    if dispersion_parameter(m)
        newrow = [_dname(m), Ryu.writefixed(dispersion(m), digits), "", "", ""]
        for rr in fnames(m)
            push!(newrow, "")
        end
        push!(rows, newrow)
    end

    tbl = Markdown.Table(rows, align)
    return tbl
end

function _markdown(s::OptSummary)
    rows = [
        ["", ""],
        ["**Initialization**", ""],
        ["Initial parameter vector", string(s.initial)],
        ["Initial objective value", string(s.finitial)],
        ["**Optimizer settings** ", ""],
        ["Optimizer (from NLopt)", "`$(s.optimizer)`"],
        ["Lower bounds", string(s.lowerbd)],
        ["`ftol_rel`", string(s.ftol_rel)],
        ["`ftol_abs`", string(s.ftol_abs)],
        ["`xtol_rel`", string(s.xtol_rel)],
        ["`xtol_abs`", string(s.xtol_abs)],
        ["`initial_step`", string(s.initial_step)],
        ["`maxfeval`", string(s.maxfeval)],
        ["`maxtime`", string(s.maxtime)],
        ["**Result**", ""],
        ["Function evaluations", string(s.feval)],
        ["Final parameter vector", "$(round.(s.final; digits=4))"],
        ["Final objective value", "$(round.(s.fmin; digits=4))"],
        ["Return code", "`$(s.returnvalue)`"],
    ]

    tbl = Markdown.Table(rows, [:l, :l])
    return tbl
end

function _markdown(vc::VarCorr)
    σρ = vc.σρ
    nmvec = string.([keys(σρ)...])
    cnmvec = string.(foldl(vcat, [keys(sig)...] for sig in getproperty.(values(σρ), :σ)))
    σvec = vcat(collect.(values.(getproperty.(values(σρ), :σ)))...)
    if !isnothing(vc.s)
        push!(σvec, vc.s)
        push!(nmvec, "Residual")
    end
    nρ = maximum(length.(getproperty.(values(σρ), :ρ)))
    varvec = abs2.(σvec)
    digits = _printdigits(σvec)
    showσvec = aligncompact(σvec, digits)
    showvarvec = aligncompact(varvec, digits)

    newrow = [" ", "Column", " Variance", "Std.Dev"]
    iszero(nρ) || push!(newrow, "Corr.")
    rows = [newrow]

    align = [:l, :l, :r, :r]
    iszero(nρ) || push!(align, :r)

    ind = 1
    for (i, v) in enumerate(values(vc.σρ))
        newrow = [string(nmvec[i])]
        firstrow = true
        k = length(v.σ)   # number of columns in grp factor k
        ρ = v.ρ
        ρind = 0
        for j in 1:k
            !firstrow && push!(newrow, " ")
            push!(newrow, string(cnmvec[ind]))
            push!(newrow, string(showvarvec[ind]))
            push!(newrow, string(showσvec[ind]))
            for l in 1:(j - 1)
                ρind += 1
                ρval = ρ[ρind]
                if ρval === -0.0
                    push!(newrow, ".")
                else
                    push!(newrow, Ryu.writefixed(ρval, 2, true))
                end
            end
            push!(rows, newrow)
            newrow = Vector{String}()
            firstrow = false
            ind += 1
        end
    end
    if !isnothing(vc.s)
        newrow = [string(last(nmvec)), " ", string(showvarvec[ind]), string(showσvec[ind])]
        push!(rows, newrow)
    end

    # pad out the rows to all have the same length
    rowlen = maximum(length, rows)
    for rr in rows
        append!(rr, repeat([" "], rowlen - length(rr)))
    end
    append!(align, repeat([:r], rowlen - length(align)))

    tbl = Markdown.Table(rows, align)
    return tbl
end


text:
"""Returns the constant forcing matrix Fx that varies only meriodionally
as a cosine with strongest forcing in the middle and vanishing forcing at the northern
and southern boundary."""
function channel_wind()
    # for non-dimensional gradients the wind forcing needs to contain the grid spacing Δ
    xx_u,yy_u = meshgrid(x_u,y_u)
    Fx = (Δ*Fx0/ρ/water_depth)*cos.(π*(yy_u/Ly .- 1/2)).^2
    return Numtype.(Fx)
end

"""Returns the constant forcing matrix Fx that varies only meriodionally
as a hyperbolic tangent with strongest shear in the middle."""
function shear_wind()
    # for non-dimensional gradients the wind forcing needs to contain the grid spacing Δ
    xx_u,yy_u = meshgrid(x_u,y_u)
    Fx = (Δ*Fx0/ρ/water_depth)*tanh.(2π*(yy_u/Ly .- 1/2))
    return Numtype.(Fx)
end

"""Returns the constant forcing matrix Fx that varies only meriodionally
with a superposition of sin & cos for a double gyre circulation.
See Cooper&Zanna 2015 or Kloewer et al 2018."""
function double_gyre_wind()
    # for non-dimensional gradients the wind forcing needs to contain the grid spacing Δ
    xx_u,yy_u = meshgrid(x_u,y_u)
    Fx = (Δ*Fx0/ρ/water_depth)*(cos.(2π*(yy_u/Ly .- 1/2)) + 2*sin.(π*(yy_u/Ly .- 1/2)))
    return Numtype.(Fx)
end

"""Returns a zero matrix for no wind forcing."""
function no_wind()
    return zeros(Numtype,nux,nuy)
end

"""Returns a reference state for Newtonian cooling/surface relaxation shaped as a
hyperbolic tangent to force the continuity equation."""
function interface_relaxation()
    # width of a tangent is defined as the distance between 65% of its minimum value and 65% of the max
    xx_T,yy_T = meshgrid(x_T,y_T)
    η_ref = -(η_refh/2)*tanh.(2π*(Ly/(4*η_refw))*(yy_T/Ly .- 1/2))
    return Numtype.(η_ref)
end

if wind_forcing == "channel"
    wind = channel_wind
elseif wind_forcing == "shear"
    wind = shear_wind
elseif wind_forcing == "double_gyre"
    wind = double_gyre_wind
elseif wind_forcing == "none"
    wind = no_wind
else
    throw(error("Wind forcing not correctly specified."))
end

max_stars_repo_name: eschnett/juls

max_issues_repo_name: eschnett/juls

max_forks_repo_name: eschnett/juls


text:
module Gameplay

using ..Wikipedia, ..Wikipedia.Articles

export newgame

const DIFFICULTY_EASY = 2
const DIFFICULTY_MEDIUM = 4
const DIFFICULTY_HARD = 6

const MAX_NUMBER_OF_STEPS = 10

function newgame(difficulty = DIFFICULTY_HARD)
  articles = Article[]

  for i in 1:difficulty+1
    article = if i == 1
                article = persistedarticle(fetchrandom()...)
              else
                url = rand(articles[i-1].links)
                existing_articles = Articles.find(url)

                article = isempty(existing_articles) ? persistedarticle(fetchpage(url)...) : existing_articles[1]
              end

    push!(articles, article)
  end

  articles
end

end

max_stars_repo_name: PacktPublishing/Julia-Projects

max_issues_repo_name: PacktPublishing/Julia-Projects

max_forks_repo_name: PacktPublishing/Julia-Projects


text:
#= Local convention for currency symbol. =#

# default: :after
const LOCAL_SYMBOL_LOCATION = Dict{Symbol, Symbol}(
    :ANG => :before,
    :AOA => :before,
    :ARS => :before,
    :AUD => :before,
    :AWG => :before,
    :CAD => :before,
    :USD => :before,
    :GBP => :before)


text:
Base.@kwdef mutable struct RateInfo
    enable::Bool = true
    num_batch::Int = 200
    t_start::UInt64 = 0
    cnt::Int = 0
    name::String = ""
end

function init(rate::RateInfo, name::String)
    rate.name = name
    if rate.cnt == 0
        rate.t_start = time_ns()
    end
end

function printrate(rate::RateInfo)
    if rate.enable && rate.cnt == rate.num_batch
        t_stop = time_ns()
        ns_elapsed = time_ns() - rate.t_start
        avg_rate = rate.num_batch / (ns_elapsed * 1e-9)
        println("Average rate of $(rate.name): ", avg_rate, " Hz")
        rate.t_start = t_stop
        rate.cnt = 0
    end
    rate.cnt += 1
end


text:
module VaxData

export AbstractVax, VaxInt, VaxFloat

abstract type AbstractVax <: Real end
abstract type VaxInt <: AbstractVax end
abstract type VaxFloat <: AbstractVax end

import Base: IEEEFloat

include("constants.jl")
include("vaxints.jl")
include("vaxfloatf.jl")
include("vaxfloatd.jl")
include("vaxfloatg.jl")
include("promote.jl")

const VaxTypes = Union{VaxInt16,VaxInt32,VaxFloatF,VaxFloatD,VaxFloatG}

function Base.read(s::IO, ::Type{T}) where T <: VaxTypes
    return read!(s, Ref{T}(0))[]::T
end

function Base.promote(x::T, y::T) where T <: AbstractVax
    Base.@_inline_meta
    px, py = Base._promote(x, y)
    Base.not_sametype((x,y), (px,py))
    px, py
end
function Base.promote(x::T, y::T, z::T) where T <: AbstractVax
    Base.@_inline_meta
    px, py, pz = Base._promote(x, y, z)
    Base.not_sametype((x,y,z), (px,py,pz))
    px, py, pz
end
function Base.promote(x::T, y::T, z::T, a::T...) where T <: AbstractVax
    p = Base._promote(x, y, z, a...)
    Base.not_sametype((x, y, z, a...), p)
    p
end

# Define common arithmetic operations (default for two of the same unknown number type is to no-op error)
# Promotion rules are such that the promotion will always be to a valid IEEE number type,
# even in the case of two identical AbstractVax types
for op in [:+, :-, :*, :/, :^, :<, :<=]
    @eval(begin
        Base.$op(x::T, y::T) where {T<:AbstractVax} = ($op)(promote(x,y)...)
    end)
end

end # module


text:
module SparseMatrixAssemblers

using SparseArrays

export SparseMatrixAssembler

"""
    SparseMatrixAssembler{Tv,Ti<:Integer}

Convenience type for assembling sparse matrices.

# Examples
```
julia> A = SparseMatrixAssembler();

julia> A[2,2] = π;

julia> Matrix(A,3,3)
3×3 Array{Float64,2}:
 0.0  0.0      0.0
 0.0  3.14159  0.0
 0.0  0.0      0.0
```
"""
struct SparseMatrixAssembler{Tv,Ti<:Integer}
    i::Vector{Ti}
    j::Vector{Ti}
    v::Vector{Tv}
end

SparseMatrixAssembler() where {Tv} = SparseMatrixAssembler{Float64}()
SparseMatrixAssembler{Tv}() where {Tv} = SparseMatrixAssembler{Tv,Int}()
SparseMatrixAssembler{Tv,Ti}() where {Tv,Ti<:Integer} = SparseMatrixAssembler{Tv,Ti}(Vector{Ti}(),Vector{Ti}(),Vector{Tv}())

Base.convert(::Type{SparseMatrixAssembler{Tv,Ti}}, A::SparseMatrixAssembler) where {Tv,Ti<:Integer} =
    SparseMatrixAssembler(
        convert(Vector{Ti},A.i),
        convert(Vector{Ti},A.j),
        convert(Vector{Tv},A.v)
    )
Base.convert(::Type{SparseMatrixAssembler{Tv}}, A::SparseMatrixAssembler) where {Tv} = convert(SparseMatrixAssembler{Tv,Int},A)

# Overload repeat to handle scalars
repeat(v::AbstractArray; kwargs...) = Base.repeat(v; kwargs...)
repeat(v::Number; kwargs...) = Base.repeat([v]; kwargs...)

function Base.setindex!(A::SparseMatrixAssembler, v, i,j)
    @assert length(v) == length(i)*length(j)
    append!(A.i, repeat(i,outer=length(j)))
    append!(A.j, repeat(j,inner=length(i)))
    append!(A.v, v)
    return v
end

SparseArrays.sparse(A::SparseMatrixAssembler) = sparse(A.i,A.j,A.v)
SparseArrays.sparse(A::SparseMatrixAssembler, m::Integer,n::Integer, combine...) = sparse(A.i,A.j,A.v,m,n, combine...)
# The following can be ambiguous with `sparse(i,j,v,m,n)`, hence we need to explicitly write out both version as above.
#     SparseArrays.sparse(A::SparseMatrixAssembler, args...) = sparse(A.i,A.j,A.v, args...)

SparseArrays.SparseMatrixCSC(A::SparseMatrixAssembler, args...) = sparse(A, args...)
Base.Matrix(A::SparseMatrixAssembler, args...) = Matrix(SparseMatrixCSC(A, args...))

end # module


text:
module SparseRegression

using LinearAlgebra
import SweepOperator: sweep!
import LearnBase: learn!, ObsDim, value, predict, deriv!
import LearningStrategies: strategy, setup!, update!, finished, cleanup!
import StatsBase: coef, AbstractWeights, Weights

using LossFunctions, PenaltyFunctions, LearningStrategies, RecipesBase

# Reexports
for pkg in [:LossFunctions, :PenaltyFunctions, :LearningStrategies]
    eval(Expr(:toplevel, Expr(:export, setdiff(names(eval(pkg)), [pkg])...)))
end

export
    SModel, ProxGrad, Fista, AdaptiveProxGrad, GradientDescent, Sweep, LinRegCholesky,
    Weights, LineSearch,
    coef, predict, value

#-----------------------------------------------------------------------# Types
abstract type Algorithm <: LearningStrategy end
abstract type GradientAlgorithm <: Algorithm end
abstract type OneIterAlgorithm <: Algorithm end
finished(a::OneIterAlgorithm, model, i) = true

include("smodel.jl")
include("algorithms/algorithms.jl")
include("algorithms/cd.jl")

#-----------------------------------------------------------------------# Auto learn!
function learn!(o::SModel; verbose::Bool = true)
    s = strategy(o)
    verbose ? learn!(o, Verbose(s)) : learn!(o, s)
    s
end

strategy(o::SModel) = strategy(AdaptiveProxGrad(o), MaxIter(), Converged(coef))

const ScaledL2 = Union{L2DistLoss, LossFunctions.ScaledDistanceLoss{L2DistLoss}}
strategy(o::SModel{<:ScaledL2, <:Union{NoPenalty, L2Penalty}}) = Sweep(o)

#-----------------------------------------------------------------------# Fake data
function fakedata(::DistanceLoss, n, p; β = range(-1, stop=1, length=p))
    x = randn(n, p)
    y = x*β + randn(n)
    x, y, collect(β)
end

function fakedata(::MarginLoss, n, p; β = range(-1, stop=1, length=p))
    x = randn(n, p)
    y = Float64[2 * (1 / (1+exp(-η)) > rand()) - 1 for η in x * β]
    x, y, collect(β)
end

end


text:
"""
    reduce(f, t::IndexedTable; select::Selection)

Apply reducer function `f` pair-wise to the selection `select` in `t`.  The reducer `f`
can be:

1. A function
1. An OnlineStat
1. A (named) tuple of functions and/or OnlineStats
1. A (named) tuple of (selector => function) or (selector => OnlineStat) pairs

# Examples

    t = table(1:5, 6:10, names = [:t, :x])

    reduce(+, t, select = :t)
    reduce((a, b) -> (t = a.t + b.t, x = a.x + b.x), t)

    using OnlineStats
    reduce(Mean(), t, select = :t)
    reduce((Mean(), Variance()), t, select = :t)

    y = reduce((min, max), t, select=:x)
    reduce((sum = +, prod = *), t, select=:x)

    # combining reduction and selection
    reduce((xsum = :x => +, negtsum = (:t => -) => +), t)
"""
function reduce(f, t::IndexedTable; select=valuenames(t), kws...)
    if haskey(kws, :init)
        return _reduce_select_init(f, t, select, kws.data.init)
    end
    _reduce_select(f, t, select)
end

function _reduce_select(f, t::Dataset, select)
    fs, input, T = init_inputs(f, rows(t, select), reduced_type, false)
    acc = init_first(fs, input[1])
    _reduce(fs, input, acc, 2)
end

function _reduce_select_init(f, t::Dataset, select, v0)
    fs, input, T = init_inputs(f, rows(t, select), reduced_type, false)
    _reduce(fs, input, v0, 1)
end

function _reduce(fs, input, acc, start)
    @inbounds @simd for i=start:length(input)
        acc = _apply(fs, acc, input[i])
    end
    acc
end

## groupreduce

addname(v, name) = v
addname(v::Tup, name::Type{<:NamedTuple}) = v
addname(v, name::Type{<:NamedTuple}) = name((v,))

struct GroupReduce{F, S, T, P, N}
    f::F
    key::S
    data::T
    perm::P
    name::N
    n::Int

    GroupReduce(f::F, key::S, data::T, perm::P; name::N = nothing) where{F, S, T, P, N} =
        new{F, S, T, P, N}(f, key, data, perm, name, length(key))
end

Base.IteratorSize(::Type{<:GroupReduce}) = Base.SizeUnknown()

function Base.iterate(iter::GroupReduce, i1=1)
    i1 > iter.n && return nothing
    f, key, data, perm, n, name = iter.f, iter.key, iter.data, iter.perm, iter.n, iter.name
    val = init_first(f, data[perm[i1]])
    i = i1+1
    while i <= n && roweq(key, perm[i], perm[i1])
        val = _apply(f, val, data[perm[i]])
        i += 1
    end
    (key[perm[i1]] => addname(val, name)), i
end

"""
    groupreduce(f, t, by = pkeynames(t); select)

Calculate a [`reduce`](@ref) operation `f` over table `t` on groups defined by the values
in selection `by`.  The result is put in a table keyed by the unique `by` values.

# Examples

    t = table([1,1,1,2,2,2], 1:6, names=[:x, :y])
    groupreduce(+,        t, :x; select = :y)
    groupreduce((sum=+,), t, :x; select = :y)  # change output column name to :sum

    t2 = table([1,1,1,2,2,2], [1,1,2,2,3,3], 1:6, names = [:x, :y, :z])
    groupreduce(+, t2, (:x, :y), select = :z)

    # different reducers for different columns
    groupreduce((sumy = :y => +, sumz = :z => +), t2, :x)
"""
function groupreduce(f, t::Dataset, by=pkeynames(t);
                     select = t isa AbstractIndexedTable ? Not(by) : valuenames(t),
                     cache=false)

    if f isa ApplyColwise
        if !(f.functions isa Union{Function, Type})
            error("Only functions are supported in ApplyColwise for groupreduce")
        end
        return groupby(grp->colwise_group_fast(f.functions, grp), t, by; select=select)
    end

    isa(f, Pair) && (f = (f,))

    data = rows(t, select)

    by = lowerselection(t, by)

    if !isa(by, Tuple)
        by=(by,)
    end
    key  = rows(t, by)
    perm = sortpermby(t, by, cache=cache)

    fs, input, T = init_inputs(f, data, reduced_type, false)

    name = isa(t, IndexedTable) ? namedtuple(nicename(f)) : nothing
    iter = GroupReduce(fs, key, input, perm, name=name)
    convert(collectiontype(t), collect_columns(iter),
            presorted=true, copy=false)
end

colwise_group_fast(f, grp::Union{Columns, Dataset}) = map(c->reduce(f, c), columns(grp))
colwise_group_fast(f, grp::AbstractVector) = reduce(f, grp)

## GroupBy

_apply_with_key(f::Tup, data::Tup, process_data) = _apply(f, map(process_data, data))
_apply_with_key(f::Tup, data, process_data) = _apply_with_key(f, columns(data), process_data)
_apply_with_key(f, data, process_data) = _apply(f, process_data(data))

_apply_with_key(f::Tup, key, data::Tup, process_data) = _apply(f, map(t->key, data), map(process_data, data))
_apply_with_key(f::Tup, key, data, process_data) = _apply_with_key(f, key, columns(data), process_data)
_apply_with_key(f, key, data, process_data) = _apply(f, key, process_data(data))

struct GroupBy
    f
    key
    data
    perm
    usekey::Bool
    name
    n::Int

    GroupBy(f, key, data, perm; usekey = false, name = nothing) =
        new(f, key, data, perm, usekey, name, length(key))
end

Base.IteratorSize(::Type{<:GroupBy}) = Base.SizeUnknown()

function Base.iterate(iter::GroupBy, i1=1)
    i1 > iter.n && return nothing
    f, key, data, perm, usekey, n, name = iter.f, iter.key, iter.data, iter.perm, iter.usekey, iter.n, iter.name
    i = i1+1
    while i <= n && roweq(key, perm[i], perm[i1])
        i += 1
    end
    process_data = t -> view(t, perm[i1:(i-1)])
    val = usekey ? _apply_with_key(f, key[perm[i1]], data, process_data) :
                   _apply_with_key(f, data, process_data)
    (key[perm[i1]] => addname(val, name)), i
end

collectiontype(::Type{<:NDSparse}) = NDSparse
collectiontype(::Type{<:IndexedTable}) = IndexedTable
collectiontype(t::Dataset) = collectiontype(typeof(t))

"""
    groupby(f, t, by = pkeynames(t); select, flatten=false)

Apply `f` to the `select`-ed columns (see [`select`](@ref)) in groups defined by the
unique values of `by`.

If `f` returns a vector, split it into multiple columns with `flatten = true`.

# Examples

    using Statistics

    t=table([1,1,1,2,2,2], [1,1,2,2,1,1], [1,2,3,4,5,6], names=[:x,:y,:z])

    groupby(mean, t, :x, select=:z)
    groupby(identity, t, (:x, :y), select=:z)
    groupby(mean, t, (:x, :y), select=:z)

    groupby((mean, std, var), t, :y, select=:z)
    groupby((q25=z->quantile(z, 0.25), q50=median, q75=z->quantile(z, 0.75)), t, :y, select=:z)

    # apply different aggregation functions to different columns
    groupby((ymean = :y => mean, zmean = :z => mean), t, :x)
"""
function groupby end

function groupby(f, t::Dataset, by=pkeynames(t);
            select = t isa AbstractIndexedTable ? Not(by) : valuenames(t),
            flatten=false, usekey = false)

    isa(f, Pair) && (f = (f,))
    data = rows(t, select)
    f = init_func(f, data)
    by = lowerselection(t, by)
    if !(by isa Tuple)
        by = (by,)
    end

    key = by == () ? fill((), length(t)) : rows(t, by)

    fs, input, S = init_inputs(f, data, reduced_type, true)

    if by == ()
        res = usekey ? _apply_with_key(fs, (), input, identity) : _apply_with_key(fs, input, identity)
        res_tup = addname(res, namedtuple(nicename(f)))
        return flatten ? res_tup[end] : res_tup
    end

    perm = sortpermby(t, by)
    # Note: we're not using S here, we'll let _groupby figure it out
    name = isa(t, IndexedTable) ? namedtuple(nicename(f)) : nothing
    iter = GroupBy(fs, key, input, perm, usekey = usekey, name = name)

    t = convert(collectiontype(t), collect_columns(iter), presorted=true, copy=false)
    t isa IndexedTable && flatten ?
        IndexedTables.flatten(t, length(columns(t))) : t
end

struct ApplyColwise{T}
    functions::T
    names
    stack::Bool
    variable::Symbol
end

ApplyColwise(f; stack = false, variable = :variable) = ApplyColwise(f, [nicename(f)], stack, variable)
ApplyColwise(t::Tuple; stack = false, variable = :variable) = ApplyColwise(t, [map(nicename,t)...], stack, variable)
ApplyColwise(t::NamedTuple; stack = false, variable = :variable) = ApplyColwise(Tuple(values(t)), keys(t), stack, variable)

init_func(f, t) = f
init_func(ac::ApplyColwise{<:Tuple}, t::AbstractVector) =
    Tuple(Symbol(n) => f for (f, n) in zip(ac.functions, ac.names))
function init_func(ac::ApplyColwise{<:Tuple}, t::Columns)
    if ac.stack
        dd -> Columns((collect(colnames(t)), ([f(x) for x in columns(dd)] for f in ac.functions)...); names = vcat(ac.variable, ac.names))
    else
        Tuple(Symbol(s, :_, n) => s => f for s in colnames(t), (f, n) in zip(ac.functions, ac.names))
    end
end

init_func(ac::ApplyColwise, t::Columns) =
    ac.stack ? dd -> Columns((collect(colnames(t)), [ac.functions(x) for x in columns(dd)]); names = vcat(ac.variable, ac.names)) :
        Tuple(s => s => ac.functions for s in colnames(t))
init_func(ac::ApplyColwise, t::AbstractVector) = ac.functions

"""
    summarize(f, t, by = pkeynames(t); select = Not(by), stack = false, variable = :variable)

Apply summary functions column-wise to a table. Return a `NamedTuple` in the non-grouped case
and a table in the grouped case. Use `stack=true` to stack results of the same summary function
for different columns.

# Examples

    using Statistics

    t = table([1, 2, 3], [1, 1, 1], names = [:x, :y])

    summarize((mean, std), t)
    summarize((m = mean, s = std), t)
    summarize(mean, t; stack=true)
    summarize((mean, std), t; select = :y)
"""
function summarize(f, t, by = pkeynames(t); select = t isa AbstractIndexedTable ? excludecols(t, by) : valuenames(t), stack = false, variable = :variable)
    flatten = stack && !(select isa Union{Int, Symbol})
    s = groupby(ApplyColwise(f, stack = stack, variable = variable), t, by, select = select, flatten = flatten)
    s isa Columns ? table(s, copy = false, presorted = true) : s
end

"""
`convertdim(x::NDSparse, d::DimName, xlate; agg::Function, vecagg::Function, name)`

Apply function or dictionary `xlate` to each index in the specified dimension.
If the mapping is many-to-one, `agg` or `vecagg` is used to aggregate the results.
If `agg` is passed, it is used as a 2-argument reduction function over the data.
If `vecagg` is passed, it is used as a vector-to-scalar function to aggregate
the data.
`name` optionally specifies a new name for the translated dimension.
"""
function convertdim(x::NDSparse, d::DimName, xlat; agg=nothing, vecagg=nothing, name=nothing, select=valuenames(x))
    ks = setcol(pkeys(x), d, d=>xlat)
    if name !== nothing
        ks = renamecol(ks, d, name)
    end

    if vecagg !== nothing
        y = convert(NDSparse, ks, rows(x, select))
        return groupby(vecagg, y)
    end

    if agg !== nothing
        return convert(NDSparse, ks, rows(x, select), agg=agg)
    end
    convert(NDSparse, ks, rows(x, select))
end

convertdim(x::NDSparse, d::Int, xlat::Dict; agg=nothing, vecagg=nothing, name=nothing, select=valuenames(x)) = convertdim(x, d, i->xlat[i], agg=agg, vecagg=vecagg, name=name, select=select)

convertdim(x::NDSparse, d::Int, xlat, agg) = convertdim(x, d, xlat, agg=agg)

sum(x::NDSparse) = sum(x.data)

"""
    reduce(f, x::NDSparse; dims)

Drop the `dims` dimension(s) and aggregate values with `f`.

    x = ndsparse((x=[1,1,1,2,2,2],
                  y=[1,2,2,1,2,2],
                  z=[1,1,2,1,1,2]), [1,2,3,4,5,6])

    reduce(+, x; dims=1)
    reduce(+, x; dims=(1,3))
"""
function Base.reduce(f, x::NDSparse; kws...)
    if haskey(kws, :dims)
        if haskey(kws, :select) || haskey(kws, :init)
            throw(ArgumentError("select and init keyword arguments cannot be used with dims"))
        end
        dims = kws.data.dims
        if dims isa Symbol
            dims = [dims]
        end
        keep = setdiff([1:ndims(x);], map(d->fieldindex(columns(x.index),d), dims))
        if isempty(keep)
            throw(ArgumentError("to remove all dimensions, use `reduce(f, A)`"))
        end
        return groupreduce(f, x, (keep...,))
    else
        select = get(kws, :select, valuenames(x))
        if haskey(kws, :init)
            return _reduce_select_init(f, x, select, kws.data.init)
        end
        return _reduce_select(f, x, select)
    end
end

"""
`reducedim_vec(f::Function, arr::NDSparse, dims)`

Like `reduce`, except uses a function mapping a vector of values to a scalar instead
of a 2-argument scalar function.
"""
function reducedim_vec(f, x::NDSparse, dims; with=valuenames(x))
    keep = setdiff([1:ndims(x);], map(d->fieldindex(columns(x.index),d), dims))
    if isempty(keep)
        throw(ArgumentError("to remove all dimensions, use `reduce(f, A)`"))
    end
    idxs, d = collect_columns(GroupBy(f, keys(x, (keep...,)), rows(x, with), sortpermby(x, (keep...,)))) |> columns
    NDSparse(idxs, d, presorted=true, copy=false)
end

reducedim_vec(f, x::NDSparse, dims::Symbol) = reducedim_vec(f, x, [dims])


text:
using CFTime
if VERSION >= v"0.7.0-beta.0"
    using Test
    using Dates
    using Printf
else
    using Base.Test
    using Compat
end

@testset "Time and calendars" begin
    include("test_time.jl")
end


text:
# include -
include("Include.jl")

max_stars_repo_name: varnerlab/BST-S-Sytem-Model-Coagulation

max_issues_repo_name: varnerlab/BST-S-Sytem-Model-Coagulation

max_forks_repo_name: varnerlab/BST-S-Sytem-Model-Coagulation


text:
@testset "MortalityTables from File" begin

    @testset "TableMetaData" begin
        @testset "whitespace management" begin
            # "2017 Loaded CSO Smoker Distinct Nonsmoker Female ALB" comes with a trailing whitespace
            @test ~isnothing(MortalityTables.table_source_map["2017 Loaded CSO Smoker Distinct Nonsmoker Female ALB"])
        end
    end

    @testset "SOA tables" begin

        @testset "Ultimate Only SOA table" begin
            cso1980 = MortalityTables.table("1980 CSO Basic Table – Male, ANB")
            @test cso1980.ultimate[35] ≈ 0.00118
            @test cso1980[35] ≈ 0.00118
            @test cso1980[35:36] ≈ [0.00118, 0.00128]
            @test cso1980[35 + 60] ≈ 0.27302
            @test cso1980[95] ≈ 0.27302
            @test_throws BoundsError cso1980[125]
            @test_throws ArgumentError survival(cso1980,10,15)
            @test_throws ArgumentError decrement(cso1980,10,15)
            @test omega(cso1980.ultimate) == 100
            @test MortalityTables.ω(cso1980) == 100
        end

        @testset "2001 VBT" begin
            vbt2001 =
                MortalityTables.table("2001 VBT Residual Standard Select and Ultimate - Male Nonsmoker, ANB")

            @test vbt2001.select[35][35] ≈ 0.00036
            @test vbt2001.select[35][35 + 60] ≈ 0.24298
            @test vbt2001.select[95][95] ≈ 0.23815
            @test vbt2001.ultimate[95] ≈ 0.24298
            @test vbt2001.ultimate[120] ≈ 1.0
            @test vbt2001.select[85][120] ≈ 1.0
            @test_throws BoundsError vbt2001.select[150]
            @test omega(vbt2001.select[20]) == 120
            @test omega(vbt2001.ultimate) == 120

            # this tests when there are ultimate rates available
            # but no select rates for the issue age
            @test_throws BoundsError vbt2001.select[116][116] ≈ 0.79988

        end

        @testset "2001 CSO" begin
            cso2001 =
                MortalityTables.table("2001 CSO Super Preferred Select and Ultimate - Male Nonsmoker, ANB")

            @test cso2001.select[35][35] ≈ 0.00037
            @test cso2001.select[35][35 + 60] ≈ 0.26719
            @test cso2001.ultimate[16] ≈ 0.00041
            @test cso2001.ultimate[95] ≈ 0.26719
            @test cso2001.ultimate[120] ≈ 1.0
            @test cso2001.select[85][120] ≈ 1.0
            @test ismissing(cso2001.select[15][15]) # age before table defines rates
            @test_throws BoundsError cso2001.select[150]
            @test_throws BoundsError cso2001.select[70][150]
            @test omega(cso2001.select[16]) == 120
            @test omega(cso2001.select[98]) == 120 # issue #71
            @test omega(cso2001.ultimate) == 120

            @testset "duration ranges" begin
                # select only
                @test all(cso2001.select[29][29:31] .== [0.00029, 0.00035, 0.0004])
                @test all(cso2001.select[35][35:59] .== [
                    0.00037,
                    0.00043,
                    0.00049,
                    0.00057,
                    0.00063,
                    0.0007,
                    0.00077,
                    0.00084,
                    0.00092,
                    0.00101,
                    0.00114,
                    0.00127,
                    0.00143,
                    0.00159,
                    0.00174,
                    0.00188,
                    0.00208,
                    0.00231,
                    0.00251,
                    0.00279,
                    0.00315,
                    0.00357,
                    0.00406,
                    0.00461,
                    0.00508,
                ])

                # crosses into ultimate
                @test all(cso2001.select[35][55:64] .== [
                    0.00315,
                    0.00357,
                    0.00406,
                    0.00461,
                    0.00508,
                    0.00621,
                    0.0069,
                    0.00773,
                    0.00867,
                    0.00965,
                ])
                @test all(cso2001.select[35][35:64] .== [
                    0.00037,
                    0.00043,
                    0.00049,
                    0.00057,
                    0.00063,
                    0.0007,
                    0.00077,
                    0.00084,
                    0.00092,
                    0.00101,
                    0.00114,
                    0.00127,
                    0.00143,
                    0.00159,
                    0.00174,
                    0.00188,
                    0.00208,
                    0.00231,
                    0.00251,
                    0.00279,
                    0.00315,
                    0.00357,
                    0.00406,
                    0.00461,
                    0.00508,
                    0.00621,
                    0.0069,
                    0.00773,
                    0.00867,
                    0.00965,
                ])
            end

        end

    end
end


text:
"""
Library of thermal system components to model heat transfer.
"""
module Thermal
using ModelingToolkit, Symbolics, IfElse, OrdinaryDiffEq

@parameters t
D = Differential(t)

include("utils.jl")

# Library of 1-dimensional heat transfer with lumped elements
include("HeatTransfer/ideal_components.jl")
include("HeatTransfer/sensors.jl")
include("HeatTransfer/sources.jl")

# Simple components for 1-dimensional incompressible thermo-fluid flow models
# TODO:
# - FluidHeatFlow

export # Interface
       HeatPort,
       # Thermal Components
       BodyRadiation, ConvectiveConductor, ConvectiveResistor,
       HeatCapacitor, ThermalConductor, ThermalResistor, ThermalCollector,
       # Thermal Sensors
       RelativeTemperatureSensor, HeatFlowSensor, TemperatureSensor,
       # Thermal Sources
       FixedHeatFlow, FixedTemperature, ThermalGround

end


text:
using Downloads: download
using OrdinaryDiffEq
using Trixi

###############################################################################
# semidiscretization of the compressible Euler equations

equations = CompressibleEulerEquations3D(1.4)

initial_condition = initial_condition_constant

boundary_conditions = Dict(
  :all => BoundaryConditionDirichlet(initial_condition)
)

solver = DGSEM(polydeg=3, surface_flux=flux_lax_friedrichs,
               volume_integral=VolumeIntegralWeakForm())

# Mapping as described in https://arxiv.org/abs/2012.12040 but reduced to 2D.
# This particular mesh is unstructured in the yz-plane, but extruded in x-direction.
# Apply the warping mapping in the yz-plane to get a curved 2D mesh that is extruded
# in x-direction to ensure free stream preservation on a non-conforming mesh.
# See https://doi.org/10.1007/s10915-018-00897-9, Section 6.
function mapping(xi, eta_, zeta_)
  # Transform input variables between -1 and 1 onto [0,3]
  eta = 1.5 * eta_ + 1.5
  zeta = 1.5 * zeta_ + 1.5

  z = zeta + 1/6 * (cos(1.5 * pi * (2 * eta - 3)/3) *
                   cos(0.5 * pi * (2 * zeta - 3)/3))

  y = eta + 1/6 * (cos(0.5 * pi * (2 * eta - 3)/3) *
                  cos(2 * pi * (2 * z - 3)/3))

  return SVector(xi, y, z)
end

# Unstructured mesh with 48 cells of the cube domain [-1, 1]^3
mesh_file = joinpath(@__DIR__, "cube_unstructured_2.inp")
isfile(mesh_file) || download("https://gist.githubusercontent.com/efaulhaber/b8df0033798e4926dec515fc045e8c2c/raw/b9254cde1d1fb64b6acc8416bc5ccdd77a240227/cube_unstructured_2.inp",
                              mesh_file)

mesh = P4estMesh{3}(mesh_file, polydeg=3,
                    mapping=mapping,
                    initial_refinement_level=0)

# Refine quadrants in y-direction of each tree at one edge to level 2
function refine_fn(p8est, which_tree, quadrant)
  if convert(Int, which_tree) < 4 && quadrant.x == 0 && quadrant.y == 0 && quadrant.level < 2
    # return true (refine)
    return Cint(1)
  else
    # return false (don't refine)
    return Cint(0)
  end
end

# Refine recursively until each desired quadrant has level 2.
# The mesh will be rebalanced before the simulation starts.
refine_fn_c = @cfunction(refine_fn, Cint, (Ptr{Trixi.p8est_t}, Ptr{Trixi.p4est_topidx_t}, Ptr{Trixi.p8est_quadrant_t}))
Trixi.refine_p4est!(mesh.p4est, true, refine_fn_c, C_NULL)

semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver, boundary_conditions=boundary_conditions)

###############################################################################
# ODE solvers, callbacks etc.

tspan = (0.0, 1.0)
ode = semidiscretize(semi, tspan)

summary_callback = SummaryCallback()

analysis_interval = 100
analysis_callback = AnalysisCallback(semi, interval=analysis_interval)

alive_callback = AliveCallback(analysis_interval=analysis_interval)

save_restart = SaveRestartCallback(interval=100,
                                   save_final_restart=true)

save_solution = SaveSolutionCallback(interval=100,
                                     save_initial_solution=true,
                                     save_final_solution=true,
                                     solution_variables=cons2prim)

stepsize_callback = StepsizeCallback(cfl=1.2)

callbacks = CallbackSet(summary_callback,
                        analysis_callback, alive_callback,
                        save_restart, save_solution,
                        stepsize_callback)

###############################################################################
# run the simulation

sol = solve(ode, CarpenterKennedy2N54(williamson_condition=false), #maxiters=1,
            dt=1.0, # solve needs some value here but it will be overwritten by the stepsize_callback
            save_everystep=false, callback=callbacks);
summary_callback() # print the timer summary