text:
###########################################################################################

# Regular grid for 1D array
function simps(y::AbstractVector{<:Number},dx::Real=1;even="avg")

	if iseven(size(y,1))
		if !(even ∈ ("avg","first","last"))
			error("even has to be one of avg, first or last."*
				"Specified even is $(even)")
		end

		T = promote_type(eltype(y), Float64)
		result = zero(T)
		val = zero(T)
		first_ind = first(axes(y,1))
		last_ind = last(axes(y,1))

		if even ∈ ["avg","first"]
			result += simpsRegular1D(view(y,first_ind:last_ind-1),dx)
			val += trapzRegular1D(view(y,last_ind-1:last_ind),dx)
		end
		if even ∈ ["avg","last"]
			result += simpsRegular1D(view(y,first_ind+1:last_ind),dx)
			val += trapzRegular1D(view(y,first_ind:first_ind+1),dx)
		end
		if even == "avg"
			val /= 2
			result /=2
		end
		result += val
	else
		result = simpsRegular1D(y,dx)
	end
	return result
end

# Regular grid for nD array
function simps(y::AbstractArray{<:Number},dx::Real=1;even="avg",axis::Integer=1)

	T = promote_type(eltype(y), Float64)

	if axis != 1
		perm = [axis]
		append!(perm,[i for i in 1:axis-1])
		append!(perm,[i for i in axis+1:ndims(y)])
		y = permutedims(y,perm)
	end

	trailing_axes,inds = trailing_indices(y)
	result = zeros(T,inds)

	if iseven(size(y,1))
		if !(even ∈ ("avg","first","last"))
			error("even has to be one of avg, first or last."*
				"Specified even is $(even)")
		end

		val = zeros(T,inds)
		first_ind = first(axes(y,1))
		last_ind = last(axes(y,1))

		if even ∈ ["avg","first"]
			simpsRegular1DAdd!(view(y,first_ind:last_ind-1,trailing_axes),
						result,dx)
			trapzRegular1DAdd!(view(y,last_ind-1:last_ind,trailing_axes),
						val,dx)
		end
		if even ∈ ["avg","last"]
			simpsRegular1DAdd!(view(y,first_ind+1:last_ind,trailing_axes),
						result,dx)
			trapzRegular1DAdd!(view(y,first_ind:first_ind+1,trailing_axes),
						val,dx)
		end
		if even == "avg"
			val ./= 2
			result ./=2
		end
		result .+= val
	else
		simpsRegular1D!(y,result,dx)
	end
	return result
end

# Irregular grid for 1D array
function simps(y::AbstractVector{<:Number},x::AbstractVector{<:Real};even="avg")

	if iseven(size(y,1))
		if !(even ∈ ("avg","first","last"))
			error("even has to be one of avg, first or last."*
				"Specified even is $(even)")
		end

		T = promote_type(eltype(y), Float64)
		result = zero(T)
		val = zero(T)
		first_ind = first(axes(y,1))
		last_ind = last(axes(y,1))
		first_ind_x = first(axes(x,1))
		last_ind_x = last(axes(x,1))

		if even ∈ ["avg","first"]
			result += simpsIrregular1D(view(y,first_ind:last_ind-1),
				view(x,first_ind_x:last_ind_x-1))
			val += trapzRegular1D(view(y,last_ind-1:last_ind),
				x[last_ind_x]-x[last_ind_x-1])
		end
		if even ∈ ["avg","last"]
			result += simpsIrregular1D(view(y,first_ind+1:last_ind),
				view(x,first_ind_x+1:last_ind_x))
			val += trapzRegular1D(view(y,first_ind:first_ind+1),
				x[first_ind_x+1]-x[first_ind_x])
		end
		if even == "avg"
			val /= 2.
			result /=2.
		end
		result += val
	else
		result = simpsIrregular1D(y,x)
	end
	return result
end

# Irregular grid for nD array
function simps(y::AbstractArray{<:Number},x::AbstractVector{<:Real};even="avg",axis::Integer=1)

	T = promote_type(eltype(y), Float64)

	if axis != 1
		perm = [axis]
		append!(perm,[i for i in 1:axis-1])
		append!(perm,[i for i in axis+1:ndims(y)])
		y = permutedims(y,perm)
	end

	trailing_axes,inds = trailing_indices(y)
	result = zeros(T,inds)

	if iseven(size(y,1))
		if !(even ∈ ("avg","first","last"))
			error("even has to be one of avg, first or last."*
				"Specified even is $(even)")
		end

		val = zeros(T,inds)

		first_ind = first(axes(y,1))
		last_ind = last(axes(y,1))

		first_ind_x = first(axes(x,1))
		last_ind_x = last(axes(x,1))

		if even ∈ ["avg","first"]
			simpsIrregular1DAdd!(view(y,first_ind:last_ind-1,trailing_axes),
				result,view(x,first_ind_x:last_ind_x-1))

			trapzRegular1DAdd!(view(y,last_ind-1:last_ind,trailing_axes),
				val,x[last_ind_x]-x[last_ind_x-1])
		end
		if even ∈ ["avg","last"]
			simpsIrregular1DAdd!(view(y,first_ind+1:last_ind,trailing_axes),
				result,view(x,first_ind_x+1:last_ind_x))

			trapzRegular1DAdd!(view(y,first_ind:first_ind+1,trailing_axes),
				val,x[first_ind_x+1]-x[first_ind_x])
		end
		if even == "avg"
			val ./= 2.
			result ./=2.
		end
		result .+= val
	else
		simpsIrregular1D!(y,result,x)
	end
	return result
end

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

function simpsweight(Nelems,index)
	@assert(1<=index<=Nelems,"index not in range")
	if index==1 || index==Nelems
		return 1
	elseif iseven(index)
		return 4
	end
	return 2
end

function quadratic_fit_integrate!(coeffs,y,x)
	# Find the lagrange polynomial fit to y(x) (quadratic)
	# integrate this between first(x) and last(x)
	@assert(length(y)==length(x)==length(coeffs)==3)

	coeffs[1] = -(x[1]-x[3])*(2x[1]-3x[2]+x[3])/(6*(x[1]-x[2]))

	coeffs[2] = -(x[1]-x[3])^3/(6*(x[1]-x[2])*(x[2]-x[3]))

	coeffs[3] = (x[1]-x[3])*(x[1]-3x[2]+2x[3])/(6*(x[2]-x[3]))

	dot(coeffs,y)
end

#############################################################################
## Simpson's 1/3 rule
#############################################################################

# 1D array with uniform grid spacing for an odd number of points
function simpsRegular1D(y::AbstractVector{<:Number},dx::Real=1)

	@assert(isodd(length(y)),"Number of elements must be odd to apply Simpson's rule")

	N = size(y,1)
	sum(simpsweight(N,ind)*yi for (ind,yi) in enumerate(y)) *dx/3
end

function simpsRegular1D!(y::AbstractVector{<:Number},
	int_y::AbstractArray{T},dx::Real=1) where {T<:Number}

	@assert(isodd(length(y)),"Number of elements must be odd to apply Simpson's rule")

	N = size(y,1)
	int_y .= sum(simpsweight(N,ind)*yi for (ind,yi) in enumerate(y)) *dx/3

	return first(int_y)
end

# UnitRanges can be integrated exactly
function simpsRegular1D(y::AbstractUnitRange{<:Real},dx::Real=1)
	(last(y)^2 - first(y)^2) * dx/2
end

# N dimensional array, integrate along the first axis
function simpsRegular1D(y::AbstractArray{<:Number},dx::Real=1)

	T = promote_type(eltype(y),Float64)
	trailing_axes = trailing_indices(y)
	int_y = zeros(T,inds)

	for ind_t in trailing_axes
		int_y[ind_t] = simpsRegular1D(view(y,:,ind_t),dx)
	end
	return int_y
end

function simpsRegular1D!(y::AbstractArray{<:Number},
	int_y::AbstractArray{<:Number},dx::Real=1)

	trailing_axes,_ = trailing_indices(y)

	for ind_t in trailing_axes
		int_y[ind_t] = simpsRegular1D(view(y,:,ind_t),dx)
	end

	return int_y
end

function simpsRegular1DAdd!(y::AbstractArray{<:Number},
	int_y::AbstractArray{<:Number},dx::Real=1)

	trailing_axes,_ = trailing_indices(y)

	for ind_t in trailing_axes
		int_y[ind_t] += simpsRegular1D(view(y,:,ind_t),dx)
	end

	return int_y
end

#############################################################################
## Simpson's 1/3 rule for irregular grids
#############################################################################

# 1D array with non-uniform grid spacing, integral along one axis
function simpsIrregular1D(y::AbstractVector{TY},x::AbstractVector{TX}) where {TY,TX}

	@assert(isodd(length(y)),"The array should have an odd number of elements")
	@assert(length(x)==length(y),"y and x need to have the same number of elements")

	T = promote_type(TY, TX)
	int_y = zero(T)

	y_no_offset = OffsetArrays.no_offset_view(y)
	x_no_offset = OffsetArrays.no_offset_view(x)

	# The strategy is to fit a quadratic to three points, and compute its integral
	coeffs = zeros(T,3)

	for ind in 1:2:length(x)-2

		int_y += quadratic_fit_integrate!(coeffs,view(y_no_offset,ind:ind+2),
						view(x_no_offset,ind:ind+2))
	end

	return int_y
end

function simpsIrregular1D!(y::AbstractVector,
	int_y::AbstractArray{<:Any,0},x::AbstractVector)

	int_y[] = simpsIrregular1D(y,x)

	return first(int_y)
end

# N dimensional array, integrate along the first axis
function simpsIrregular1D(y::AbstractArray{TY},x::AbstractVector{TX}) where {TY,TX}

	T = promote_type(TY,TX)
	trailing_axes,inds = trailing_indices(y)
	int_y = zeros(T,inds)

	for ind_t in trailing_axes
		int_y[ind_t] = simpsIrregular1D(view(y,:,ind_t),x)
	end

	return int_y
end

function simpsIrregular1D!(y::AbstractArray{<:Number},
	int_y::AbstractArray{T},x::AbstractVector{<:Real}) where {T<:Number}

	trailing_axes,_ = trailing_indices(y)

	for ind_t in trailing_axes
		int_y[ind_t] = simpsIrregular1D(view(y,:,ind_t),x)
	end

	return int_y
end

function simpsIrregular1DAdd!(y::AbstractArray{<:Number},
	int_y::AbstractArray{T},x::AbstractVector{<:Real}) where {T<:Number}

	trailing_axes,_ = trailing_indices(y)

	for ind_t in trailing_axes
		int_y[ind_t] += simpsIrregular1D(view(y,:,ind_t),x)
	end

	return int_y
end

function simpsCircle(y::AbstractVector{<:Number},
	ϕ::AbstractArray{<:Real,1}=LinRange(0,2π,size(y,1));
	even::String="avg",ϕ2π::Bool=true)::Array{<:Number,0}

	if ϕ2π
		z = simps(y,x=ϕ,even=even)
	else

		info("Assuming that 2π is not included,"*
			" will append it to the array before integrating.\n"*
			"If 2π is included set ϕ2π=true in the function call")

		check_last_value(ϕ,2π) # check if the last value of ϕ is 2π

		push!(y,y[1])

		z = simps(y,x=push!(ϕ,2π),even=even)
	end

	return z
end

function simpsSpherical2D(y::AbstractArray{<:Number},
	θ::AbstractArray{<:Real,1}=LinRange(0,π,size(y,1)),
	ϕ::AbstractArray{<:Real,1}=LinRange(0,2π,size(y,2));
	even="avg",θint="clenshaw_quadrature",GLweights=nothing,
	ϕ2π=true)::Array{<:Number,0}

	if θint=="gauss_quadrature"

		if GLweights != nothing
			z = gauss_quadrature(y,GLweights)
		end

	elseif θint=="clenshaw_quadrature"

		z = clenshaw_curtis_quadrature(y)

	else
		z = simps(sin.(θ).*y,x=θ,even=even)
	end

	simpsCircle(z,ϕ,even=even,ϕ2π=ϕ2π)
end

function simpsSpherical3D(y::AbstractArray{<:Number},
	r::Union{Real,AbstractArray{<:Real,1}}=1,
	θ::AbstractArray{<:Real,1}=LinRange(0,π,size(y,2)),
	ϕ::AbstractArray{<:Real,1}=LinRange(0,π,size(y,3));
	even="avg",θint="clenshaw_quadrature",GLweights = nothing,
	ϕ2π=true)

	if typeof(r) <: AbstractArray
		z = simps(r.^2 .* y,x=r,even=even)
	else
		dr = r
		z = dr^2*simps(@. ((1:size(y,1))^2 * y),dx=dr,even=even)
	end

	simpsSpherical2D(z,θ,ϕ,θint=θint,
		GLweights=GLweights,even=even,ϕ2π=ϕ2π)
end


text:
"""
    WindowedAssociativeOp{T,Op,Op!,V<:AbstractVector{T}}

State associated with a windowed aggregation of a binary associative operator.

If `Op!` is not specifiied, it will default to `Op`.  However, for non-bitstypes, it can be
beneficial to provide this methods to reduce memory allocations.
`V` will default to a `Vector{T}`.

# Method

Wherever summation is discussed, we can consider any alternative binary, associative,
operator. For example: `+, *, max, min, &&, union`

NB. It is interesting to observe that commutativity is *not* required by this algorithm.

Conceptually the window is maintained in two buffers:

        [---- A ---)[----- B ------)
            <                      >    <-- current window finishes at the end of B, and
                                            starts somewhere in A.

`A` is stored as a sequence of cumulative sums, such that as the "<" advances we merely pick
out the correct element:

        x_i,   x_i-1 + x_i,  x_i-2 + x_i-1 + x_i

`B` is stored as both:
- The sequence of values seen:  `x_i+1,  x_i+2,  x_i+3,  ...`
- The total of that sequence:  `x_i+1 + x_i+2 + x_i+3 + ...`

When the "<" advances from `A` to `B`, we discard `A`, and the subset of `B` remaining after
`<` becomes the new `A`. In becoming `A`, we transform its representation into that of the
cumulative sums. We create a new, empty, `B`.

`O(1)` amortized runtime complexity, and `O(L)` space complexity, where `L` is the typical
window length.

# Type parameters
- `T`: The type of the values of the array.
- `Op`: Any binary, associative, function.
- `Op!`: Op!(x, y) will perform `x + y`, storing the result in `x`.
- `V`: The subtype of AbstractVector{T} used for internal state.

# Fields
- `previous_cumsum::Vector{T}`: Corresponds to array `A` above.
- `ri_previous_cumsum::Int`: A reverse index into `previous_cumsum`, once it contains
    values. It should be subtracted from `end` in order to obtain the appropriate index.
- `values::Vector{T}`: Corresponds to array `B` above.
- `sum::T`: The sum of the elements in values.
"""
mutable struct WindowedAssociativeOp{T,Op,Op!,V<:AbstractVector{T}}
    previous_cumsum::V
    ri_previous_cumsum::Int
    values::V
    sum::T  # Will start uninitialised.

    function WindowedAssociativeOp{T,Op,Op!,V}(
        previous_cumsum::V, values::V
    ) where {T,Op,Op!,V<:AbstractVector{T}}
        return new{T,Op,Op!,V}(previous_cumsum, 0, values)
    end
end

function WindowedAssociativeOp{T,Op,Op!}() where {T,Op,Op!}
    return WindowedAssociativeOp{T,Op,Op!,Vector{T}}(T[], T[])
end
WindowedAssociativeOp{T,Op}() where {T,Op} = WindowedAssociativeOp{T,Op,Op}()

"""
    update_state!(
        state::WindowedAssociativeOp,
        value,
        num_dropped_from_window::Integer
    ) -> state

Add the specified value to the state, drop some number of elements from the start of the
window, and return `state` (which will have been mutated).

# Arguments
- `state::WindowedAssociativeOp`: The state to update (will be mutated).
- `value`: The value to add to the end of the window - must be convertible to a `T`.
- `num_dropped_from_window::Integer`: The number of elements to remove from the front of
    the window.

# Returns
- The instance `state` that was passed in.
"""
Base.@propagate_inbounds function update_state!(
    state::WindowedAssociativeOp{T,Op,Op!}, value, num_dropped_from_window::Integer
) where {T,Op,Op!}
    # Our index into previous_cumsum is advanced by the number of values we drop from the
    # window.
    state.ri_previous_cumsum += num_dropped_from_window

    # If this has taken us out-of-range of the _previous_cumsum values, we must recompute
    # them.
    elements_remaining = length(state.previous_cumsum) - state.ri_previous_cumsum

    if elements_remaining < 0
        # We have overshot the end of previous_cumsum.

        # We may also need to discard elements from values. This could happen if
        # num_dropped_from_window > 1.
        num_values_to_remove = -elements_remaining
        @boundscheck if num_values_to_remove > length(state.values)
            throw(
                ArgumentError(
                    "num_dropped_from_window = $num_dropped_from_window is out of range"
                ),
            )
        end

        # We now generate the partial sum, and set our index back to zero. values is also
        # emptied, since its information is now reflected in previous_cumsum.
        # NOTE: We need to take care here in the case of non-commutation. In accumulate, we
        # will be getting:
        #
        #    (x0, op(x0, x1), op(op(x0, x1), x2), ...)
        #
        # but we actually want:
        #
        #    (x0, op(x1, x0), op(x2, op(x1, x0)), ...)

        # Conceptually the following code is equivalent to the following, however avoids
        # unnecessary allocations:
        # trimmed_reversed_values = state.values[end:-1:1 + num_values_to_remove]
        # state.previous_cumsum = accumulate(
        #     (x, y) -> Op(y, x),
        #     trimmed_reversed_values
        # )

        empty!(state.previous_cumsum)
        upper = length(state.values)
        lower = 1 + num_values_to_remove
        if (upper - lower) >= 0  # i.e. we have a non-zero range
            i = upper
            accumulation = @inbounds(state.values[i])
            while true
                push!(state.previous_cumsum, accumulation)
                i -= 1
                i >= lower || break
                # If we were to use a mutating operation here, then we'd have to introduce
                # a copy above. So there is no drawback to using the non-mutating Op.
                accumulation = Op(@inbounds(state.values[i]), accumulation)
            end
        end

        state.ri_previous_cumsum = 0
        empty!(state.values)
        # state.sum is now garbage, but we are not going to use it before we recompute it.
    end

    # We want to make sure that we never actually mutate any value object that we pass to
    # `update_state!`.
    # Note that when initialising `state.sum`, we should take a copy if our mutating `Op!`
    # is different to `Op`. This is because we know that `Op` is necessarily non-mutating.
    _copy(x) = (Op == Op!) ? x : deepcopy(x)

    # Include the new value in sum and values.
    state.sum = isempty(state.values) ? _copy(value) : Op!(state.sum, value)
    push!(state.values, value)
    return state
end

"""
    window_value(state::WindowedAssociativeOp{T})::T where T

Get the value currently represented by the state.

Behaviour is undefined if this is called when the window is empty.

# Arguments:
- `state::WindowedAssociativeOp{T}`: The state to query.

# Returns:
- `T`: The result of aggregating over the values in the window.
"""
function window_value(state::WindowedAssociativeOp{T,Op})::T where {T,Op}
    # Include contributions both from A and B buffers.
    # Remember that we are indexing from the back.
    index = length(state.previous_cumsum) - state.ri_previous_cumsum
    return if index == 0
        # We aren't using the A buffer, either because values is full or the A buffer has
        # not yet been populated.
        state.sum
    else
        Op(@inbounds(state.previous_cumsum[index]), state.sum)
    end
end

"""
    function window_size(state::WindowedAssociativeOp)::Int

Get the current size of the window in `state`.

# Arguments:
- `state::WindowedAssociativeOp`: The state to query.

# Returns:
- `Int`: The current size of the window.
"""
function window_size(state::WindowedAssociativeOp)::Int
    # Include contributions both from A and B buffers.
    # Remember that we are indexing from the back.
    index = length(state.previous_cumsum) - state.ri_previous_cumsum
    return length(state.values) + index
end

ext: jl


text:
using PiGPIO

red_pin1 = 18 #change these numbers accordingly with your GPIO pins
red_pin2 = 23

p=Pi()

set_mode(p, red_pin1, PiGPIO.OUTPUT)
set_mode(p, red_pin2, PiGPIO.OUTPUT)

try
    for i in 1:10
        PiGPIO.write(p, red_pin1, PiGPIO.HIGH)
        PiGPIO.write(p, red_pin2, PiGPIO.LOW)
        sleep(0.5)
        PiGPIO.write(p, red_pin1, PiGPIO.LOW)
        PiGPIO.write(p, red_pin2, PiGPIO.HIGH)
        sleep(0.5)
    end
finally
    println("Cleaning up!")
    set_mode(p, red_pin1, PiGPIO.INPUT)
    set_mode(p, red_pin2, PiGPIO.INPUT)
end


text:
@inline function flow(::Type{Val{:exproy}}, θ::AbstractVector{T}, σ::T,    logp::T, ψ::T, d::Integer, d1::Integer, Dgt0::Bool, sgn_ext::Bool, geoid::Real, roy::T) where {T}
    if d == 0
        sgn_ext && return θ[9]
        return zero(T)
    end
    u = rev_exp(θ[1],θ[2],θ[3],θ[4],θ[5],σ,logp,ψ,Dgt0,geoid, roy) + (d==1 ?  θ[6] : θ[7] ) # + θ[8]*d)
    d>1      && (u *= d)
    d1 == 1  && (u += θ[8])
    return u::T
end

@inline function flow(::Type{Val{:exproy_Dgt0}}, θ::AbstractVector{T}, σ::T,    logp::T, ψ::T, d::Integer, d1::Integer, Dgt0::Bool, sgn_ext::Bool, geoid::Real, roy::T) where {T}
    if d == 0
        sgn_ext && return θ[9]
        return zero(T)
    end
    u = rev_exp(θ[1],θ[2],θ[3],θ[4],θ[5],σ,logp,ψ,Dgt0,geoid, roy) + (!Dgt0 ? θ[6] : θ[7] ) # + θ[8]*d)
    d>1      && (u *= d)
    d1 == 1  && (u += θ[8])
    return u::T
end

@inline function flow(::Type{Val{:exproy_extend}}, θ::AbstractVector{T}, σ::T,    logp::T, ψ::T, d::Integer, d1::Integer, Dgt0::Bool, sgn_ext::Bool, geoid::Real, roy::T) where {T}
    if d == 0
        sgn_ext && return θ[9] + θ[10]*exp(ψ)
        return zero(T)
    end
    u = rev_exp(θ[1],θ[2],θ[3],θ[4],θ[5],σ,logp,ψ,Dgt0,geoid, roy) + (d==1 ?  θ[6] : θ[7] ) # + θ[8]*d)
    d>1      && (u *= d)
    d1 == 1  && (u += θ[8])
    return u::T
end

@inline function flow(::Type{Val{:exp}}, θ::AbstractVector{T}, σ::T,    logp::T, ψ::T, d::Integer, d1::Integer, Dgt0::Bool, sgn_ext::Bool, geoid::Real, roy::T) where {T}
    if d == 0
        sgn_ext && return θ[8]
        return zero(T)
    end
    u = rev_exp(1,θ[1],θ[2],θ[3],θ[4],σ,logp,ψ,Dgt0,geoid, roy) + (d==1 ?  θ[5] : θ[6] ) # + θ[8]*d)
    d>1      && (u *= d)
    d1 == 1  && (u += θ[7])
    return u::T
end

@inline function flow(::Type{Val{:d1_priordrill}}, θ::AbstractVector{T}, σ::T,    logp::T, ψ::T, d::Integer, d1::Integer, Dgt0::Bool, sgn_ext::Bool, geoid::Real, roy::T) where {T}
    if d == 0
        sgn_ext && return θ[7]
        return zero(T)
    end
    u = rev_exp(1,θ[1],1,θ[2],θ[3],σ,logp,ψ,Dgt0,geoid, roy) + (d==1 ?  θ[4] : θ[5] ) # + θ[8]*d)
    d>1      && (u *= d)
    d1 == 1  && (u += θ[6])
    return u::T
end

@inline function flow(::Type{Val{:Dgt0_priordrill}}, θ::AbstractVector{T}, σ::T,    logp::T, ψ::T, d::Integer, d1::Integer, Dgt0::Bool, sgn_ext::Bool, geoid::Real, roy::T) where {T}
    if d == 0
        sgn_ext && return θ[7]
        return zero(T)
    end
    u = rev_exp(1,θ[1],1,θ[2],θ[3],σ,logp,ψ,Dgt0,geoid, roy) + (d==1 ?  θ[4] : θ[5] ) # + θ[8]*d)
    d>1      && (u *= d)
    d1 == 1  && (u += θ[6])
    return u::T
end

@inline function flow(::Type{Val{:d1cost_extend}}, θ::AbstractVector{T}, σ::T,    logp::T, ψ::T, d::Integer, d1::Integer, Dgt0::Bool, sgn_ext::Bool, geoid::Real, roy::T) where {T}
    if d == 0
        sgn_ext && return θ[7] + θ[8]*exp(ψ)
        return zero(T)
    end
    u = rev_exp(1,θ[1],1,θ[2],θ[3],σ,logp,ψ,Dgt0,geoid, roy) + (d==1 ?  θ[4] : θ[5] ) # + θ[8]*d)
    d>1      && (u *= d)
    d1 == 1  && (u += θ[6])
    return u::T
end

@inline function flow(::Type{Val{:onecost_extend}}, θ::AbstractVector{T}, σ::T,    logp::T, ψ::T, d::Integer, d1::Integer, Dgt0::Bool, sgn_ext::Bool, geoid::Real, roy::T) where {T}
    if d == 0
        sgn_ext && return θ[5] + θ[6]*exp(ψ)
        return zero(T)
    end
    u = rev_exp(1,θ[1],1,θ[2],θ[3],σ,logp,ψ,Dgt0,geoid, roy) + θ[4]
    d>1      && (u *= d)
    return u::T
end

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

@inline function flowdθ(::Type{Val{:exproy}}, θ::AbstractVector{T}, σ::T,     logp::T, ψ::T, k::Integer,d::Integer,           d1::Integer, Dgt0::Bool, sgn_ext::Bool, geoid::Real, roy::T)::T where {T}
    d == 0 && !sgn_ext && return zero(T)

    # revenue
    k == 1  && return - d * exp(θ[2] + θ[3]*logp + θ[4]*geoid + Eexpψ(θ[5], σ, ψ, Dgt0) ) * roy
    k == 2  && return   d * rev_exp(θ[1],θ[2],θ[3],θ[4],θ[5],σ,logp,ψ,Dgt0,geoid, roy)
    k == 3  && return   d * rev_exp(θ[1],θ[2],θ[3],θ[4],θ[5],σ,logp,ψ,Dgt0,geoid, roy) * logp
    k == 4  && return   d * rev_exp(θ[1],θ[2],θ[3],θ[4],θ[5],σ,logp,ψ,Dgt0,geoid, roy) * convert(T,geoid)
    k == 5  && return   d * rev_exp(θ[1],θ[2],θ[3],θ[4],θ[5],σ,logp,ψ,Dgt0,geoid, roy) * ( Dgt0 ? ψ : ψ*_ρ(σ) + θ[k]*(1-_ρ2(σ)))

    # drilling cost
    k == 6  && return  d  == 1 ? one(T)  : zero(T)
    k == 7  && return  d  == 1 ? zero(T) : convert(T,d)
    # k == 8  && return  d  == 1 ? zero(T) : convert(T,d^2)
    k == 8  && return  d1 == 1 ? one(T)  : zero(T)

    # extension cost
    k == 9  && return d == 0 && sgn_ext ? one(T) : zero(T)
    # k == 10 && return d == 0 && sgn_ext ? exp(ψ) : zero(T)

    throw(error("$k out of bounds"))
end

@inline function flowdθ(::Type{Val{:exproy_Dgt0}}, θ::AbstractVector{T}, σ::T,     logp::T, ψ::T, k::Integer,d::Integer,           d1::Integer, Dgt0::Bool, sgn_ext::Bool, geoid::Real, roy::T)::T where {T}
    d == 0 && !sgn_ext && return zero(T)

    # revenue
    k == 1  && return - d * exp(θ[2] + θ[3]*logp + θ[4]*geoid + Eexpψ(θ[5], σ, ψ, Dgt0) ) * roy
    k == 2  && return   d * rev_exp(θ[1],θ[2],θ[3],θ[4],θ[5],σ,logp,ψ,Dgt0,geoid, roy)
    k == 3  && return   d * rev_exp(θ[1],θ[2],θ[3],θ[4],θ[5],σ,logp,ψ,Dgt0,geoid, roy) * logp
    k == 4  && return   d * rev_exp(θ[1],θ[2],θ[3],θ[4],θ[5],σ,logp,ψ,Dgt0,geoid, roy) * convert(T,geoid)
    k == 5  && return   d * rev_exp(θ[1],θ[2],θ[3],θ[4],θ[5],σ,logp,ψ,Dgt0,geoid, roy) * ( Dgt0 ? ψ : ψ*_ρ(σ) + θ[k]*(1-_ρ2(σ)))

    # drilling cost
    k == 6  && return  !Dgt0   ? convert(T,d) : zero(T)
    k == 7  && return  !Dgt0   ? zero(T)      : convert(T,d)
    k == 8  && return  d1 == 1 ? one(T)       : zero(T)

    # extension cost
    k == 9  && return d == 0 && sgn_ext ? one(T) : zero(T)

    throw(error("$k out of bounds"))
end

@inline function flowdθ(::Type{Val{:exproy_extend}}, θ::AbstractVector{T}, σ::T,     logp::T, ψ::T, k::Integer,d::Integer,           d1::Integer, Dgt0::Bool, sgn_ext::Bool, geoid::Real, roy::T)::T where {T}
    d == 0 && !sgn_ext && return zero(T)

    # revenue
    k == 1  && return - d * exp(θ[2] + θ[3]*logp + θ[4]*geoid + Eexpψ(θ[5], σ, ψ, Dgt0) ) * roy
    k == 2  && return   d * rev_exp(θ[1],θ[2],θ[3],θ[4],θ[5],σ,logp,ψ,Dgt0,geoid, roy)
    k == 3  && return   d * rev_exp(θ[1],θ[2],θ[3],θ[4],θ[5],σ,logp,ψ,Dgt0,geoid, roy) * logp
    k == 4  && return   d * rev_exp(θ[1],θ[2],θ[3],θ[4],θ[5],σ,logp,ψ,Dgt0,geoid, roy) * convert(T,geoid)
    k == 5  && return   d * rev_exp(θ[1],θ[2],θ[3],θ[4],θ[5],σ,logp,ψ,Dgt0,geoid, roy) * ( Dgt0 ? ψ : ψ*_ρ(σ) + θ[k]*(1-_ρ2(σ)))

    # drilling cost
    k == 6  && return  d  == 1 ? one(T)  : zero(T)
    k == 7  && return  d  == 1 ? zero(T) : convert(T,d)
    k == 8  && return  d1 == 1 ? one(T)  : zero(T)

    # extension cost
    k == 9  && return d == 0 && sgn_ext ? one(T) : zero(T)
    k == 10 && return d == 0 && sgn_ext ? exp(ψ) : zero(T)

    throw(error("$k out of bounds"))
end

@inline function flowdθ(::Type{Val{:exp}}, θ::AbstractVector{T}, σ::T,     logp::T, ψ::T, k::Integer,d::Integer,           d1::Integer, Dgt0::Bool, sgn_ext::Bool, geoid::Real, roy::T)::T where {T}
    d == 0 && !sgn_ext && return zero(T)

    # revenue
    k == 1  && return   d * rev_exp(1,θ[1],θ[2],θ[3],θ[4],σ,logp,ψ,Dgt0,geoid, roy)
    k == 2  && return   d * rev_exp(1,θ[1],θ[2],θ[3],θ[4],σ,logp,ψ,Dgt0,geoid, roy) * logp
    k == 3  && return   d * rev_exp(1,θ[1],θ[2],θ[3],θ[4],σ,logp,ψ,Dgt0,geoid, roy) * convert(T,geoid)
    k == 4  && return   d * rev_exp(1,θ[1],θ[2],θ[3],θ[4],σ,logp,ψ,Dgt0,geoid, roy) * ( Dgt0 ? ψ : ψ*_ρ(σ) + θ[k]*(1-_ρ2(σ)))

    # drilling cost
    k == 5  && return  d  == 1 ? one(T)  : zero(T)
    k == 6  && return  d  == 1 ? zero(T) : convert(T,d)
    k == 7  && return  d1 == 1 ? one(T)  : zero(T)

    # extension cost
    k == 8  && return d == 0 && sgn_ext ? one(T) : zero(T)

    throw(error("$k out of bounds"))
end

@inline function flowdθ(::Type{Val{:extend_constr_onecost}}, θ::AbstractVector{T}, σ::T,     logp::T, ψ::T, k::Integer,d::Integer,           d1::Integer, Dgt0::Bool, sgn_ext::Bool, geoid::Real, roy::T)::T where {T}
    d == 0 && !sgn_ext && return zero(T)

    # revenue
    k == 1  && return   d * rev_exp(1,θ[1],1,θ[2],θ[3],σ,logp,ψ,Dgt0,geoid,roy)
    k == 2  && return   d * rev_exp(1,θ[1],1,θ[2],θ[3],σ,logp,ψ,Dgt0,geoid,roy) * geoid
    k == 3  && return   d * rev_exp(1,θ[1],1,θ[2],θ[3],σ,logp,ψ,Dgt0,geoid,roy) * ( Dgt0 ? ψ : ψ*_ρ(σ) + θ[k]*(1-_ρ2(σ)))

    # drilling cost
    k == 4  && return d == 0 ? zero(T) : convert(T,d)

    # extension cost
    k == 5  && return d == 0 && sgn_ext ? one(T) : zero(T)
    k == 6  && return d == 0 && sgn_ext ? exp(ψ) : zero(T)

    throw(error("$k out of bounds"))
end

@inline function flowdθ(::Type{Val{:exproy_extend_constr}}, θ::AbstractVector{T}, σ::T,     logp::T, ψ::T, k::Integer,d::Integer,           d1::Integer, Dgt0::Bool, sgn_ext::Bool, geoid::Real, roy::T)::T where {T}
    d == 0 && !sgn_ext && return zero(T)

    # revenue
    k == 1  && return   d * rev_exp(1,θ[1],1,θ[2],θ[3],σ,logp,ψ,Dgt0,geoid, roy)
    k == 2  && return   d * rev_exp(1,θ[1],1,θ[2],θ[3],σ,logp,ψ,Dgt0,geoid, roy) * convert(T,geoid)
    k == 3  && return   d * rev_exp(1,θ[1],1,θ[2],θ[3],σ,logp,ψ,Dgt0,geoid, roy) * ( Dgt0 ? ψ : ψ*_ρ(σ) + θ[k]*(1-_ρ2(σ)))

    # drilling cost
    k == 4  && return  d  == 1 ? one(T)  : zero(T)
    k == 5  && return  d  == 1 ? zero(T) : convert(T,d)
    k == 6  && return  d1 == 1 ? one(T)  : zero(T)

    # extension cost
    k == 7  && return d == 0 && sgn_ext ? one(T) : zero(T)
    k == 8  && return d == 0 && sgn_ext ? exp(ψ) : zero(T)

    throw(error("$k out of bounds"))
end

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

@inline function flowdσ(::FF, θ::AbstractVector{T}, σ::T, logp::T, ψ::T, d::Integer, geoid::Real, roy::T)::T where {FF <: Union{Type{Val{:exproy}}, Type{Val{:exproy_extend}}, Type{Val{:exproy_Dgt0}}}, T}
    if d == 0
        return zero(T)
    else
        return d * drevdσ_exp(θ[1],θ[2],θ[3],θ[4],θ[5],σ,logp,ψ,geoid, roy)
    end
end

@inline function flowdσ(::Type{Val{:exp}}, θ::AbstractVector{T}, σ::T, logp::T, ψ::T, d::Integer, geoid::Real, roy::T)::T where {T}
    if d == 0
        return zero(T)
    else
        return d * drevdσ_exp(1,θ[1],θ[2],θ[3],θ[4],σ,logp,ψ,geoid, roy)
    end
end

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

@inline function flowdψ(::FF, θ::AbstractVector{T}, σ::T, logp::T, ψ::T, d::Integer, sgn_ext::Bool, geoid::Real, roy::T) where {T, FF <: Union{Type{Val{:exproy}}, Type{Val{:exproy_Dgt0}}}}
    if d == 0
        return zero(T) # sgn_ext ? θ[10] : zero(T)
    else
        return (d * drevdψ_exp(θ[1],θ[2],θ[3],θ[4],θ[5],σ,logp,ψ,geoid, roy))::T
    end
end

@inline function flowdψ(::Type{Val{:exproy_extend}}, θ::AbstractVector{T}, σ::T, logp::T, ψ::T, d::Integer, sgn_ext::Bool, geoid::Real, roy::T)::T where {T}
    if d == 0
        return sgn_ext ? θ[10]*exp(ψ) : zero(T)
    else
        return d * drevdψ_exp(θ[1],θ[2],θ[3],θ[4],θ[5],σ,logp,ψ,geoid, roy)
    end
end

@inline function flowdψ(::Type{Val{:exp}}, θ::AbstractVector{T}, σ::T, logp::T, ψ::T, d::Integer, sgn_ext::Bool, geoid::Real, roy::T) where {T}
    if d == 0
        return zero(T) # sgn_ext ? θ[10] : zero(T)
    else
        return (d * drevdψ_exp(1,θ[1],θ[2],θ[3],θ[4],σ,logp,ψ,geoid, roy))::T
    end
end

@inline function flowdψ(::Type{Val{:exproy_extend_constr}}, θ::AbstractVector{T}, σ::T, logp::T, ψ::T, d::Integer, sgn_ext::Bool, geoid::Real, roy::T)::T where {T}
    if d == 0
        return sgn_ext ? θ[8]*exp(ψ) : zero(T)
    else
        return d * drevdψ_exp(1,θ[1],1,θ[2],θ[3],σ,logp,ψ,geoid, roy)
    end
end

@inline function flowdψ(::Type{Val{:extend_constr_onecost}}, θ::AbstractVector{T}, σ::T, logp::T, ψ::T, d::Integer, sgn_ext::Bool, geoid::Real, roy::T)::T where {T}
    if d == 0
        return sgn_ext ? θ[6]*exp(ψ) : zero(T)
    else
        return d * drevdψ_exp(1,θ[1],1,θ[2],θ[3],σ,logp,ψ,geoid, roy)
    end
end


text:
# import of Rmath functions

module RFunctions

import Rmath: libRmath

### import macro

function _import_rmath(rname::Symbol, jname::Symbol, pargs)
    # C function names
    if rname == :norm
        dfun = Expr(:quote, "dnorm4")
        pfun = Expr(:quote, "pnorm5")
        qfun = Expr(:quote, "qnorm5")
        rfun = Expr(:quote, "rnorm")
    else
        dfun = Expr(:quote, string('d', rname))    # density
        pfun = Expr(:quote, string('p', rname))    # cumulative probability
        qfun = Expr(:quote, string('q', rname))    # quantile
        rfun = Expr(:quote, string('r', rname))    # random sampling
    end

    # Julia function names
    pdf = Symbol(jname, "pdf")
    cdf = Symbol(jname, "cdf")
    ccdf = Symbol(jname, "ccdf")

    logpdf = Symbol(jname, "logpdf")
    logcdf = Symbol(jname, "logcdf")
    logccdf = Symbol(jname, "logccdf")

    invcdf = Symbol(jname, "invcdf")
    invccdf = Symbol(jname, "invccdf")
    invlogcdf = Symbol(jname, "invlogcdf")
    invlogccdf = Symbol(jname, "invlogccdf")

    is_tukey = rname == :tukey

    rand = Symbol(jname, "rand")
    has_rand = true
    if rname == :nbeta || rname == :nf || rname == :nt || is_tukey
        has_rand = false
    end

    # arguments & argument types
    _pts = fill(:Cdouble, length(pargs))

    if is_tukey
        dtypes = Expr(:tuple, :Cdouble, :Cdouble, _pts..., :Cint)
        ptypes = Expr(:tuple, :Cdouble, :Cdouble, _pts..., :Cint, :Cint)
        qtypes = Expr(:tuple, :Cdouble, :Cdouble, _pts..., :Cint, :Cint)
        rtypes = Expr(:tuple, :Cdouble, _pts...)
    else
        dtypes = Expr(:tuple, :Cdouble, _pts..., :Cint)
        ptypes = Expr(:tuple, :Cdouble, _pts..., :Cint, :Cint)
        qtypes = Expr(:tuple, :Cdouble, _pts..., :Cint, :Cint)
        rtypes = Expr(:tuple, _pts...)
    end

    pdecls = [Expr(:(::), ps, :(Union{Float64,Int})) for ps in pargs] # [:(p1::Union{Float64, Int}), :(p2::Union{...}), ...]

    if is_tukey
        # ptukey and qtukey have an extra literal 1 argument
        pargs = (1, pargs...)
    end

    # Function implementation
    quote
        if $(!is_tukey)
            $pdf($(pdecls...), x::Union{Float64,Int}) =
                ccall(($dfun, libRmath), Float64, $dtypes, x, $(pargs...), 0)

            $logpdf($(pdecls...), x::Union{Float64,Int}) =
                ccall(($dfun, libRmath), Float64, $dtypes, x, $(pargs...), 1)
        end

        $cdf($(pdecls...), x::Union{Float64,Int}) =
            ccall(($pfun, libRmath), Float64, $ptypes, x, $(pargs...), 1, 0)

        $ccdf($(pdecls...), x::Union{Float64,Int}) =
            ccall(($pfun, libRmath), Float64, $ptypes, x, $(pargs...), 0, 0)

        $logcdf($(pdecls...), x::Union{Float64,Int}) =
            ccall(($pfun, libRmath), Float64, $ptypes, x, $(pargs...), 1, 1)

        $logccdf($(pdecls...), x::Union{Float64,Int}) =
            ccall(($pfun, libRmath), Float64, $ptypes, x, $(pargs...), 0, 1)

        $invcdf($(pdecls...), q::Union{Float64,Int}) =
            ccall(($qfun, libRmath), Float64, $qtypes, q, $(pargs...), 1, 0)

        $invccdf($(pdecls...), q::Union{Float64,Int}) =
            ccall(($qfun, libRmath), Float64, $qtypes, q, $(pargs...), 0, 0)

        $invlogcdf($(pdecls...), lq::Union{Float64,Int}) =
            ccall(($qfun, libRmath), Float64, $qtypes, lq, $(pargs...), 1, 1)

        $invlogccdf($(pdecls...), lq::Union{Float64,Int}) =
            ccall(($qfun, libRmath), Float64, $qtypes, lq, $(pargs...), 0, 1)

        if $has_rand
            $rand($(pdecls...)) =
                ccall(($rfun, libRmath), Float64, $rtypes, $(pargs...))
        end
    end
end

macro import_rmath(rname, jname, pargs...)
    esc(_import_rmath(rname, jname, pargs))
end

### Import specific functions

@import_rmath beta beta α β
@import_rmath binom binom n p
@import_rmath chisq chisq k
@import_rmath f fdist ν1 ν2
@import_rmath gamma gamma α β
@import_rmath hyper hyper ms mf n
@import_rmath norm norm μ σ
@import_rmath nbinom nbinom r p
@import_rmath pois pois λ
@import_rmath t tdist k
@import_rmath tukey srdist k ν

@import_rmath nbeta nbeta α β λ
@import_rmath nchisq nchisq k λ
@import_rmath nf nfdist k1 k2 λ
@import_rmath nt ntdist k λ

end


text:
module TestTwoInertiasAndIdealGear

using TinyModia
using DifferentialEquations
using ModiaPlot
using Unitful

TwoInertiasAndIdealGearTooManyInits = Model(
    J1 = 0.0025,
    J2 = 170,
    r  = 105,
    tau_max = 1,
    phi1 = Var(init = 0.0),
    w1   = Var(init = 1.0),
    phi2 = Var(init = 0.5),
    w2   = Var(init = 0.0),
    equations = :[
        tau = if time < 1u"s"; tau_max elseif time < 2u"s"; 0 elseif time < 3u"s"; -tau_max else 0 end,

        # inertia1
        w1 = der(phi1),
        J1*der(w1) = tau - tau1,

        # gear
        phi1   = r*phi2,
        r*tau1 = tau2,

        # inertia2]
        w2 = der(phi2),
        J2*der(w2) = tau2
    ]
)

TwoInertiasAndIdealGear = TwoInertiasAndIdealGearTooManyInits | Map(phi1 = Var(init=nothing), w1=Var(init=nothing))

twoInertiasAndIdealGearTooManyInits = @instantiateModel(TwoInertiasAndIdealGearTooManyInits)
twoInertiasAndIdealGear             = @instantiateModel(TwoInertiasAndIdealGear)

println("Next simulate! should result in an error:\n")
simulate!(twoInertiasAndIdealGearTooManyInits, Tsit5(), stopTime = 4.0, log=true)

simulate!(twoInertiasAndIdealGear, Tsit5(), stopTime = 4.0, log=false,
          logParameters=true, logStates=true,
          requiredFinalStates=[1.5628074713622309, -6.878080753044174e-5])

plot(twoInertiasAndIdealGear, ["phi2", "w2"])

end


text:
#SafeTestsets makes every test run in a separate enviromment
#so that errors can be found independently
using SafeTestsets
@safetestset "Operators fixed tests" begin include("operators_fixed_test.jl") end
@safetestset "States fixed tests" begin include("states_fixed_test.jl") end


text:
"""
`numparameters(f)`

Returns the number of parameters of `f` for the method which has the most parameters.
"""
function numargs(f)
  typ = Tuple{Any, Val{:analytic}, Vararg}
  typ2 = Tuple{Any, Type{Val{:analytic}}, Vararg} # This one is required for overloaded types
  numparam = maximum([(m.sig<:typ || m.sig<:typ2) ? 0 : length(m.sig.parameters) for m in methods(f)]) #in v0.5, all are generic
  return (numparam-1) #-1 in v0.5 since it adds f as the first parameter
end

function isinplace(f,inplace_param_number)
  numargs(f)>=inplace_param_number
end

function isinplace{iip}(f::AbstractParameterizedFunction{iip},inplace_param_number)
  iip
end

macro def(name, definition)
    return quote
        macro $(esc(name))()
            esc($(Expr(:quote, definition)))
        end
    end
end

macro muladd(ex)
  esc(to_muladd(ex))
end

function to_muladd(ex)
  is_add_operation(ex) || return ex

  all_operands = ex.args[2:end]
  mul_operands = filter(is_mul_operation, all_operands)
  odd_operands = filter(x->!is_mul_operation(x), all_operands)

  muladd_operands = collect(zip(
    to_muladd.((x->x.args[2]).(mul_operands)),
    to_muladd.((x->x.args[3]).(mul_operands))))

  if isempty(odd_operands)
    to_be_muladded = muladd_operands[1:end-1]
    last_operation = :($(muladd_operands[end][1]) * $(muladd_operands[end][2]))
  else
    to_be_muladded = muladd_operands
    last_operation = make_addition(odd_operands)
  end

  foldr(last_operation, to_be_muladded) do xs, r
    :($(Base.muladd)($(xs[1]), $(xs[2]), $r))
  end
end

is_operation(ex::Expr, op::Symbol) = ex.head == :call && !isempty(ex.args) && ex.args[1] == op
is_operation(ex, op::Symbol) = false

is_add_operation(ex) = is_operation(ex, :+)
is_mul_operation(ex) = is_operation(ex, :*)

make_addition(args) = length(args) == 1 ? args[1] : Expr(:call, :+, args...)


text:
module ChainRules

using Reexport
@reexport using ChainRulesCore

using Base.Broadcast: materialize, materialize!, broadcasted, Broadcasted, broadcastable
using LinearAlgebra
using LinearAlgebra.BLAS
using Random
using Requires
using Statistics

# Basically everything this package does is overloading these, so we make an exception
# to the normal rule of only overload via `ChainRulesCore.rrule`.
import ChainRulesCore: rrule, frule

if VERSION < v"1.3.0-DEV.142"
    # In prior versions, the BLAS submodule also exported `dot`, which caused a conflict
    # with its parent module. To get around this, we can simply create a hard binding for
    # the one we want to use without qualification.
    import LinearAlgebra: dot
end

include("rulesets/Core/core.jl")

include("rulesets/Base/utils.jl")
include("rulesets/Base/nondiff.jl")
include("rulesets/Base/base.jl")
include("rulesets/Base/fastmath_able.jl")
include("rulesets/Base/evalpoly.jl")
include("rulesets/Base/array.jl")
include("rulesets/Base/arraymath.jl")
include("rulesets/Base/mapreduce.jl")

include("rulesets/Statistics/statistics.jl")

include("rulesets/LinearAlgebra/utils.jl")
include("rulesets/LinearAlgebra/blas.jl")
include("rulesets/LinearAlgebra/dense.jl")
include("rulesets/LinearAlgebra/structured.jl")
include("rulesets/LinearAlgebra/factorization.jl")

include("rulesets/Random/random.jl")

# Note: The following is only required because package authors sometimes do not
# declare their own rules using `ChainRulesCore.jl`. For arguably good reasons.
# So we define them here for them.
function __init__()
    @require NaNMath="77ba4419-2d1f-58cd-9bb1-8ffee604a2e3" begin
        include("rulesets/packages/NaNMath.jl")
    end

    @require SpecialFunctions="276daf66-3868-5448-9aa4-cd146d93841b" begin
        include("rulesets/packages/SpecialFunctions.jl")
    end
end

end # module


text:
include("rosenbrock.jl")

x0 = ones(6)
c(x) = [sum(x)]
nlp2 = ADNLPModel(rosenbrock,  x0,
                 lvar = fill(-10.0,size(x0)), uvar = fill(10.0,size(x0)),
                 y0 = [0.0], c = c, lcon = [-Inf], ucon = [6.])

nlp_at_x_c = NLPAtX(x0, NaN*ones(nlp2.meta.ncon))
stop_nlp_c = NLPStopping(nlp2, KKT, nlp_at_x_c)

a = zeros(6)
fill_in!(stop_nlp_c, a)

@test cons(nlp2, a) == stop_nlp_c.current_state.cx
@test jac(nlp2, a) == stop_nlp_c.current_state.Jx

@test stop!(stop_nlp_c) == false
# we make sure the counter of stop works properly
@test stop_nlp_c.meta.nb_of_stop == 1

sol = ones(6)
fill_in!(stop_nlp_c, sol)

@test stop!(stop_nlp_c) == true
@test stop_nlp_c.current_state.current_score == KKT(stop_nlp_c.pb, stop_nlp_c.current_state)

stop_nlp_default = NLPStopping(nlp2, atol = 1.0)
fill_in!(stop_nlp_default, sol)
@test stop_nlp_default.meta.atol == 1.0
@test stop!(stop_nlp_default) == true

stop_nlp_default.optimality_check = unconstrained_check
update_and_stop!(stop_nlp_default, cx = nothing)
@test stop_nlp_default.current_state.cx != nothing

#Keywords in the stop! call
nlp_at_x_kargs = NLPAtX(x0, NaN*ones(nlp2.meta.ncon))
stop_nlp_kargs = NLPStopping(nlp2, (x,y; test = 1.0, kwargs...) -> Stopping.KKT(x,y,kwargs...) + test, nlp_at_x_c)
fill_in!(stop_nlp_kargs, sol)
@test stop!(stop_nlp_kargs) == false
@test stop!(stop_nlp_kargs, test = 0.0) == true


text:
@testset "innerjoin" begin
    l = 2:5
    r = 3:6
    @test innerjoin(l, r) == [3, 4, 5]
    @test innerjoin(l, [3,3,3]) == [3, 3, 3]
    @test innerjoin(identity, identity, l, r) == [3, 4, 5]
    @test innerjoin(identity, identity, tuple, l, r) == [(3,3), (4,4), (5,5)]
    @test innerjoin(identity, identity, tuple, isequal, l, r) == [(3,3), (4,4), (5,5)]
    @test innerjoin(identity, identity, tuple, ==, l, r) == [(3,3), (4,4), (5,5)]

    @test innerjoin(identity, identity, tuple, isless, l, r) == [(2,3), (2,4), (2,5), (2,6), (3,4), (3,5), (3,6), (4,5), (4,6), (5,6)]
end


text:
using Mycelia
using Documenter
import Weave
import DataStructures

PKG_BASE = dirname(dirname(pathof(Mycelia)))

mkpath("$(PKG_BASE)/docs/src")
# copy readme from repo to being main index file in documentation
cp("$(PKG_BASE)/README.md", "$(PKG_BASE)/docs/src/index.md", force=true)

for documentation_group in filter(x -> isdir(x) && !occursin(r"^\.", basename(x)), readdir("$(PKG_BASE)/docs/_src", join=true))
    # convert all current .ipynb notebooks to markdown files
#     @show documentation_group
#     @show basename(documentation_group)
    outdir = joinpath("$(PKG_BASE)/docs/src", basename(documentation_group))
    mkpath(outdir)
    if basename(documentation_group) == "5.Development"
        for notebook in filter(x -> occursin(r".ipynb$", x), readdir(documentation_group, join=true))
            out_path = joinpath("$(PKG_BASE)/docs/src", basename(documentation_group), basename(notebook) * ".md")
            Weave.convert_doc(notebook, out_path)
        end
    #     copy non notebook files over
        for non_notebook in filter(x -> !occursin(r".ipynb$", x), readdir(documentation_group, join=true))
            out_path = joinpath("$(PKG_BASE)/docs/src", basename(documentation_group), basename(non_notebook))
            cp(non_notebook, out_path, force=true)
        end
    else
        for notebook in filter(x -> occursin(r".ipynb$", x), readdir(documentation_group, join=true))
            out_path = joinpath("$(PKG_BASE)/docs/src", basename(documentation_group), basename(notebook) * ".md")
            Weave.weave(notebook, out_path=out_path, doctype="github")
        end
    #     copy non notebook files over
        for non_notebook in filter(x -> !occursin(r".ipynb$", x), readdir(documentation_group, join=true))
            out_path = joinpath("$(PKG_BASE)/docs/src", basename(documentation_group), basename(non_notebook))
            cp(non_notebook, out_path, force=true)
        end
    end
end

pages = DataStructures.SortedDict()
for documentation_group in filter(x -> isdir(x) && !occursin(r"^\.", basename(x)), readdir("$(PKG_BASE)/docs/src", join=true))
    md_pages = filter(x -> occursin(r"\.md", x), readdir(documentation_group))
    md_pages = map(page -> "$(basename(documentation_group))/$(page)", md_pages)
    pages[basename(documentation_group)] = md_pages
end

pages=[
    "Home" => "index.md",
    collect(pages)...
]

makedocs(;
    modules=[Mycelia],
    authors="Cameron Prybol <cameron.prybol@gmail.com> and contributors",
    repo="https://github.com/cjprybol/Mycelia/blob/{commit}{path}#L{line}",
    sitename="Mycelia",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://cjprybol.github.io/Mycelia",
        assets=String[],
        collapselevel=1,
    ),
    pages=pages
)

deploydocs(
    repo = "github.com/cjprybol/Mycelia.git",
    push_preview = true,
    deps = nothing,
    make = nothing,
    devurl = "docs"
)

# rm("$(PKG_BASE)/docs/src", recursive=true)

max_stars_repo_name: snystrom/Mycelia

max_issues_repo_name: snystrom/Mycelia

max_forks_repo_name: cjprybol/Eisenia


text:
push!(LOAD_PATH, realpath(pwd() * "/../src"))

using Aminda, ImageView

vs = view_specs(:dino_sparse_ring)
voi = bounding_box(:dino_sparse_ring)
vl = View(vs[1])
v = View(vs[2])
vr = View(vs[3])
nbrs = Set{View}([vl, vr])

dm = DepthMap(v, nbrs, voi)

view(dm.depth, pixelspacing = [1,1])

max_stars_repo_name: twadleigh/aminda

max_issues_repo_name: twadleigh/aminda

max_forks_repo_name: twadleigh/aminda


text:
using AbstractAlgebra

@testset "isgroebner certify" begin
    R, (x, y, z) = PolynomialRing(GF(2^31 - 1), ["x","y","z"], ordering=:degrevlex)

    @test Groebner.isgroebner([R(1)], certify=true)
    @test Groebner.isgroebner([x], certify=true)
    @test Groebner.isgroebner([x, y], certify=true)
    @test Groebner.isgroebner([x, x], certify=true)
    @test Groebner.isgroebner([x, y, z], certify=true)

    @test !Groebner.isgroebner([x + y, x], certify=true)
    @test Groebner.isgroebner([z*x, z*x, R(1)], certify=true)

    fs = Groebner.change_ordering(Groebner.rootn(3, ground=GF(2^31-1)), :degrevlex)
    @test !Groebner.isgroebner(fs, certify=true)
    @test Groebner.isgroebner(Groebner.groebner(fs), certify=true)

    fs = Groebner.change_ordering(Groebner.rootn(3, ground=QQ), :degrevlex)
    @test !Groebner.isgroebner(fs, certify=true)
    @test Groebner.isgroebner(Groebner.groebner(fs), certify=true)

    fs = Groebner.change_ordering(Groebner.noonn(2, ground=GF(2^31-1)), :degrevlex)
    @test !Groebner.isgroebner(fs, certify=true)
    @test Groebner.isgroebner(Groebner.groebner(fs), certify=true)

    fs = Groebner.change_ordering(Groebner.noonn(2, ground=QQ), :degrevlex)
    @test !Groebner.isgroebner(fs, certify=true)
    @test Groebner.isgroebner(Groebner.groebner(fs), certify=true)

    fs = Groebner.change_ordering(Groebner.noonn(6, ground=GF(2^31-1)), :degrevlex)
    @test !Groebner.isgroebner(fs, certify=true)
    @test Groebner.isgroebner(Groebner.groebner(fs), certify=true)

    fs = Groebner.change_ordering(Groebner.noonn(6, ground=QQ), :degrevlex)
    @test !Groebner.isgroebner(fs, certify=true)
    @test Groebner.isgroebner(Groebner.groebner(fs), certify=true)
end


text:
w = [1.0, 0.1, 0.3, 1.0]
cosa = cos(0.5)
sina = sin(0.5)
KitBase.global_frame(randn(2), cosa, sina)
KitBase.local_frame(randn(2), cosa, sina)
KitBase.global_frame(w, cosa, sina)
KitBase.local_frame(w, cosa, sina)

KitBase.global_frame(rand(3), rand(3, 3))
KitBase.local_frame(rand(3), rand(3, 3))
KitBase.global_frame(rand(5), rand(3, 3))
KitBase.local_frame(rand(5), rand(3, 3))

KitBase.unit_normal(rand(2), rand(2))
KitBase.unit_normal(rand(3), rand(3), rand(3))
KitBase.point_distance(rand(2), rand(2), rand(2))

#--- structure mesh ---#
KitBase.PSpace1D()
KitBase.PSpace1D(0.0, 1.0) |> show
KitBase.PSpace2D()
KitBase.PSpace2D(0.0, 1.0, 0.0, 1.0) |> show
KitBase.CSpace2D(0.0, 1.0, 10, 0.0, π, 10) |> show

x0 = 0.0
nx = 15
dx = 0.1
KitBase.uniform_mesh(x0, nx, dx)

x = collect(0:0.1:1)
y = collect(0:0.1:1)
z = collect(0:0.1:1)
KitBase.ndgrid(x)
KitBase.ndgrid(x, y)
KitBase.ndgrid(x, y, z)
KitBase.meshgrid(x)
KitBase.meshgrid(x, y)
KitBase.meshgrid(x, y, z)

KitBase.find_idx(randn(20), 0.13, mode = :uniform)
KitBase.find_idx(randn(20), 0.13, mode = :nonuniform)

#--- unstructure mesh ---#
cd(@__DIR__)
ps = KitBase.UnstructPSpace("t1.msh")
KitBase.write_vtk(ps.points, ps.cellid, randn(size(ps.cellid, 1)))


text:
function construct_ERK_coefficients(stages)

    if stages == 2
        alpha=zeros(stages,stages)
        beta=zeros(stages,1)
        alpha[2,1]=1
        beta[:]=alpha[stages,:]
    elseif stages == 3
        alpha=zeros(stages,stages)
        beta=zeros(stages,1)
        alpha[2,1]=0.5857864376269049511983112757903019214303281246230519
        alpha[3,1]=0.02859547920896831706610375859676730714344270820768398
        alpha[3,2]=0.9714045207910316829338962414032326928565572917923160
        beta[1]=0.3535533905932737622004221810524245196424179688442370
        beta[2]=0.3535533905932737622004221810524245196424179688442370
        beta[3]=0.2928932188134524755991556378951509607151640623115260
    elseif stages == 4
        alpha=zeros(stages,stages)
        beta=zeros(stages,1)
        alpha[2,1]=1767732205903/2027836641118
        alpha[3,1]=5535828885825/10492691773637
        alpha[3,2]=788022342437/10882634858940
        alpha[4,1]=6485989280629/16251701735622
        alpha[4,2]=-4246266847089/9704473918619
        alpha[4,3]=10755448449292/10357097424841
        beta[1]=1471266399579/7840856788654
        beta[2]=-4482444167858/7529755066697
        beta[3]=11266239266428/11593286722821
        beta[4]=1767732205903/4055673282236
    elseif stages == 6
        alpha=zeros(stages,stages)
        beta=zeros(stages,1)
        alpha[2,1]=1.0/2.0
        alpha[3,1]=13861/62500
        alpha[3,2]=6889/62500
        alpha[4,1]=-116923316275/2393684061468
        alpha[4,2]=-2731218467317/15368042101831
        alpha[4,3]=9408046702089/11113171139209
        alpha[5,1]=-451086348788/2902428689909
        alpha[5,2]=-2682348792572/7519795681897
        alpha[5,3]=12662868775082/11960479115383
        alpha[5,4]=3355817975965/11060851509271
        alpha[6,1]=647845179188/3216320057751
        alpha[6,2]=73281519250/8382639484533
        alpha[6,3]=552539513391/3454668386233
        alpha[6,4]=3354512671639/8306763924573
        alpha[6,5]=4040/17871
        beta[1]=82889/524892
        beta[2]=0
        beta[3]=15625/83664
        beta[4]=69875/102672
        beta[5]=-2260/8211
        beta[6]=1/4
    elseif stages == 8
        alpha=zeros(stages,stages)
        beta=zeros(stages,1)
        alpha[2,1]=41/100
        alpha[3,1]=367902744464/2072280473677
        alpha[3,2]=677623207551/8224143866563
        alpha[4,1]=1268023523408/10340822734521
        alpha[4,2]=0
        alpha[4,3]=1029933939417/13636558850479
        alpha[5,1]=14463281900351/6315353703477
        alpha[5,2]=0
        alpha[5,3]=66114435211212/5879490589093
        alpha[5,4]=-54053170152839/4284798021562
        alpha[6,1]=14090043504691/34967701212078
        alpha[6,2]=0
        alpha[6,3]=15191511035443/11219624916014
        alpha[6,4]=-18461159152457/12425892160975
        alpha[6,5]=-281667163811/9011619295870
        alpha[7,1]=19230459214898/13134317526959
        alpha[7,2]=0
        alpha[7,3]=21275331358303/2942455364971
        alpha[7,4]=-38145345988419/4862620318723
        alpha[7,5]=-1/8
        alpha[7,6]=-1/8
        alpha[8,1]=-19977161125411/11928030595625
        alpha[8,2]=0
        alpha[8,3]=-40795976796054/6384907823539
        alpha[8,4]=177454434618887/12078138498510
        alpha[8,5]=782672205425/8267701900261
        alpha[8,6]=-69563011059811/9646580694205
        alpha[8,7]=7356628210526/4942186776405
        beta[1]=-872700587467/9133579230613
        beta[2]=0
        beta[3]=0
        beta[4]=22348218063261/9555858737531
        beta[5]=-1143369518992/8141816002931
        beta[6]=-39379526789629/19018526304540
        beta[7]=32727382324388/42900044865799
        beta[8]=41/200
    elseif stages == 10
        alpha=zeros(stages,stages)
        beta=zeros(stages,1)
        alpha[2,1]=0.2928932188134525
        alpha[3,1]=0.3602847895715037
        alpha[3,2]=0.02099677689026724
        alpha[4,1]=0.4267095308101442
        alpha[4,3]=0.04296038329994516
        alpha[5,1]=0.4901881062634202
        alpha[5,4]=0.06787015549498755
        alpha[6,1]=0.548902422602367
        alpha[6,5]=0.09754418680435926
        alpha[7,1]=0.6003576019251088
        alpha[7,6]=0.1344773551299359
        alpha[8,1]=0.6404229388942057
        alpha[8,7]=0.1828003658091574
        alpha[9,1]=0.6612274743186617
        alpha[9,8]=0.2503841780330199
        alpha[10,1]=0.3942428457889902
        alpha[10,9]=0.6057571542110098
        beta[2]=0.7071067811865476
        beta[10]=0.2928932188134525
    end

    return alpha, beta
end


text:
using MiniBee
using Plots
using TrajectoryOptimization
using StaticArrays
using Sockets
using Altro

model = Spacecraft(kt=1, max_thrust=10000, max_amplitude_angular_control=100)

tf = 10.
N = 100
xf_pos = SA[100, 0, 100.]
waypoints = Waypoint.([
    SA[60, 0., 40],
    SA[75, 0., 80],
], [33, 66])   # in knot points
dt = repeat([tf/(N-1)], N-1)
u0 = [copy(SA[1000.0, 0, 1., 0.]) for _ ∈ 1:N-1]

opts = SolverOptions(
    penalty_scaling=5.,
    penalty_initial=10.,
    constraint_tolerance=1e-5,
    # penalty_max=
)

traj = main(model, opts, N, tf, dt, xf_pos, u0)

write(model, traj, "mini_bee_data.json")

plot(traj.time, traj.position', title="(t₀=0, tf=$tf) - $N points", legend=:topleft)
plot(traj.time[1:end-1], traj.thrust', title="Thrust")
plot(traj.time[1:end-1], traj.angular_controls', title="Angular controls")


text:
module TestDPO
using Test
using Catlab.Graphs
using Catlab.Present
using Catlab.WiringDiagrams
using Catlab.CategoricalAlgebra.FinSets
using Catlab, Catlab.Theories, Catlab.CategoricalAlgebra
using Catlab.Theories: AttrDesc

# Wiring diagrams
#################

A, B, C, D = [:A], [:B], [:C], [:D];
f,g,h,i = Box(:f, A, B), Box(:g, B, C),Box(:h, C, D), Box(:i, [:D,:B], [:D,:B]);

GWD = WiringDiagram([:A,:A],[:D,:B])
fv, gv, hv, fv2 = add_box!(GWD, f), add_box!(GWD, g), add_box!(GWD, h), add_box!(GWD, f);
add_wires!(GWD, Pair[
  (input_id(GWD),1) => (fv,1),
  (fv,1) => (gv,1),
  (gv,1) => (hv,1),
  (hv,1) => (output_id(GWD),1),
  (input_id(GWD),2) => (fv2,1),
  (fv2,1) => (output_id(GWD),2)]);

RWD = WiringDiagram([], []);
fv, iv, gv = add_box!(RWD, Box(:f,[],B)), add_box!(RWD, i), add_box!(RWD, Box(:g,B,[]));
add_wires!(RWD, Pair[(iv,1) => (iv, 1),
                     (fv, 1) => (iv, 2),
                     (iv, 2)=>(gv,1)]);

LWD = WiringDiagram([], []);
fv, gv = add_box!(LWD, Box(:f,[],B)), add_box!(LWD, Box(:g,B,[]));
add_wires!(LWD, Pair[
  (fv,1) => (gv,1),
]);

IWD = WiringDiagram([], []);
fv, gv = add_box!(IWD, Box(:f,[],[])), add_box!(IWD, Box(:g,[],[]));

XWD = WiringDiagram([:A,:A],[:D,:B])
fv, gv, hv = add_box!(XWD, f), add_box!(XWD, g), add_box!(XWD, h)
fv2, iv = add_box!(XWD, f), add_box!(XWD, i);
add_wires!(XWD, Pair[
  (input_id(XWD),1) => (fv,1),
  (iv,1) => (iv, 1), (fv, 1) => (iv, 2), (iv, 2)=>(gv,1),
  (gv,1) => (hv,1),
  (hv,1) => (output_id(XWD),1),
  (input_id(XWD),2) => (fv2,1),
  (fv2,1) => (output_id(XWD),2)]);

L=ACSetTransformation(IWD.diagram, LWD.diagram, Box=[1,2]);
R=ACSetTransformation(IWD.diagram, RWD.diagram, Box=[1,3]);
m=ACSetTransformation(LWD.diagram, GWD.diagram, Box=[1,2],InPort=[2],OutPort=[1],Wire=[1]);
@test is_isomorphic(XWD.diagram, rewrite_match(L,R,m))

# Graphs with attributes
########################

@present TheoryDecGraph(FreeSchema) begin
  E::Ob
  V::Ob
  src::Hom(E,V)
  tgt::Hom(E,V)

  X::Data
  dec::Attr(E,X)
end

@present TheoryLabelledDecGraph <: TheoryDecGraph begin
  label::Attr(V,X)
end

const LabelledDecGraph = ACSetType(TheoryLabelledDecGraph, index=[:src,:tgt])

aI2 = @acset LabelledDecGraph{String} begin
  V = 2
  E = 0
  label = ["a","b"]
end

aarr = @acset LabelledDecGraph{String} begin
  V = 2
  E = 1
  src = [1]
  tgt = [2]
  dec = ["e1"]
  label = ["a","b"]
end

abiarr = @acset LabelledDecGraph{String} begin
  V = 2
  E = 2
  src = [1,2]
  tgt = [2,1]
  dec = ["e1","rev_e1"]
  label = ["a","b"]
end

aspan = @acset LabelledDecGraph{String} begin
  V = 3
  E = 2
  src = [1,1]
  tgt = [2,3]
  dec = ["e1","e2"]
  label = ["a","b","c"]
end

expected = @acset LabelledDecGraph{String} begin
  V = 3
  E = 3
  src = [1,1,2]
  tgt = [2,3,1]
  dec = ["e1","e2","rev_e1"]
  label = ["a","b","c"]
end

L = ACSetTransformation(aI2, aarr, V=[1,2]);
R = ACSetTransformation(aI2, abiarr, V=[1,2]);
m = ACSetTransformation(aarr, aspan, V=[2,1], E=[1]);  # sends 'a'->'b' and 'b'->'a'

@test_throws ErrorException("ACSet colimit does not exist: label attributes a != b") rewrite_match(L,R,m)

m = ACSetTransformation(aarr, aspan, V=[1,2], E=[1]);

@test is_isomorphic(expected, rewrite_match(L, R, m))

# Graphs
########

# Example graphs
I2 = Graph(2)
I3 = Graph(3)
#   e1   e2
# 1 <- 2 -> 3
span = Graph(3)
add_edges!(span, [2,2],[1,3])
# 1 -> 2
arr = path_graph(Graph, 2)
# 1 <-> 2
biarr = Graph(2)
add_edges!(biarr, [1,2],[2,1])
# 1 -> 2 -> 3 -> 1
tri = cycle_graph(Graph, 3)
# 4 <- 1 -> 2 and 2 <- 3 -> 4
dispan = Graph(4)
add_edges!(dispan, [1,1,3,3],[2,4,2,4])

#      e1
#    1 -> 2
# e2 |  / ^
#    vv   | e4
#    4 -> 3
#      e5
squarediag = Graph(4)
add_edges!(squarediag, [1,1,2,3,4],[2,4,4,2,3])

# Add a reverse arrow to span
span_w_arrow = Graph(3)
add_edges!(span_w_arrow,[1,1,2],[2,3,1])

L = CSetTransformation(I2, arr, V=[1,2])
R = CSetTransformation(I2, biarr, V=[1,2])
m = CSetTransformation(arr, span, V=[2,1], E=[1])
@test is_isomorphic(span_w_arrow, rewrite_match(L, R, m))

# Remove apex of a subspan (top left corner of squarediag, leaves the triangle behind)
L = CSetTransformation(I2, span, V=[1,3])
m = CSetTransformation(span, squarediag, V=[2,1,4], E=[1,2])
@test is_isomorphic(tri, rewrite_match(L,id(I2),m))

# Remove self-edge using a *non-monic* match morphism
two_loops = Graph(2)
add_edges!(two_loops,[1,2],[1,2]) # ↻1   2↺
one_loop = Graph(2)
add_edges!(one_loop,[2],[2]) # 1   2↺

L = CSetTransformation(I2, arr, V=[1,2])
m = CSetTransformation(arr, two_loops, V=[1, 1], E=[1])
@test is_isomorphic(one_loop, rewrite_match(L,id(I2),m))

# Simplest non-trivial, non-monic exmaple
@present TheoryFinSet(FreeSchema) begin
  X::Ob
end
const FinSetType = ACSetType(TheoryFinSet)

I, L, G = [FinSetType() for _ in 1:3]
add_parts!(I,:X,2)
add_parts!(L,:X,1)
add_parts!(G,:X,1)
l = CSetTransformation(I,L,X=[1,1])
m = CSetTransformation(L,G,X=[1])
@test isempty(dangling_condition(l,m))
@test isempty(vcat(collect(id_condition(l,m))...))
ik, kg = pushout_complement(l,m)
# There are 3 functions `ik` that make this a valid P.C.
# codom=1 with [1,1], codom=2 with [1,2] or [2,1]
K = ik.codom
@test nparts(K, :X) == 1 # algorithm currently picks the first option

# Non-discrete interface graph. Non-monic matching
arr_loop= Graph(2)
add_edges!(arr_loop,[1,2],[2,2])  # 1->2↺
arrarr = Graph(2)
add_edges!(arrarr, [1,1],[2,2]) #  1⇉2
arrarr_loop = Graph(2)
add_edges!(arrarr_loop,[1,1,2],[2,2,2]) # 1⇉2↺
arr_looploop = Graph(2)
add_edges!(arr_looploop, [1,2,2],[2,2,2]) # 1-> ↻2↺

L = CSetTransformation(arr, arr, V=[1,2],E=[1]) # identity
R = CSetTransformation(arr, arrarr, V=[1,2], E=[1])
m = CSetTransformation(arr, arr_loop, V=[2,2], E=[2]) # NOT MONIC
@test is_isomorphic(arr_looploop, rewrite_match(L,R,m))

# only one monic match
@test is_isomorphic(arrarr_loop, rewrite(L, R, arr_loop, true))

# two possible morphisms L -> squarediag, but both violate dangling condition
L = CSetTransformation(arr, span, V=[1,2], E=[1]);
m = CSetTransformation(span, squarediag, V=[2,1,4], E=[1,2]);

@test (:src, 5, 4) in dangling_condition(L,m)

# violate id condition because two orphans map to same point
L = CSetTransformation(I2, biarr, V=[1,2]); # delete both arrows
m = CSetTransformation(biarr, arr_loop, V=[2,2], E=[2,2]);
@test (:E, 1, 2, 2) in id_condition(L,m)[1]
L = CSetTransformation(arr, biarr, V=[1,2], E=[1]); # delete one arrow
@test (:E, 1, 2) in id_condition(L,m)[2]

span_triangle = Graph(3); # 2 <- 1 -> 3 (with edge 2->3)
add_edges!(span_triangle,[1,1,2], [2,3,3]);

L = CSetTransformation(arr, tri, V=[1,2], E=[1]);
m = CSetTransformation(tri, squarediag, V=[2,4,3], E=[3,5,4]);
@test is_isomorphic(span_triangle, rewrite_match(L,id(arr),m))

k, g = pushout_complement(L, m); # get PO complement to do further tests

# the graph interface is equal to the final graph b/c we only delete things
@test is_isomorphic(span_triangle, codom(k))

# Check pushout properties 1: apex is the original graph
@test is_isomorphic(squarediag, pushout(L, k).cocone.apex) # recover original graph

# Check pushout properties 2: the diagram commutes
Lm = compose(L,m);
kg = compose(k,g);
for I_node in 1:2
  @test Lm[:V](I_node) == kg[:V](I_node)
end
@test Lm[:E](1) == kg[:E](1)

# Check pushout properties 3: for every pair of unmatched things between K and L, they are NOT equal
for K_node in 1:3
  @test m[:V](3) != g[:V](K_node)
end

for K_edge in 2:3
  @test m[:E](3) != g[:E](K_edge)
end

# Undirected bipartite graphs
#############################

# 1 --- 1
#    /
# 2 --- 2

z_ = UndirectedBipartiteGraph()
add_vertices₁!(z_, 2)
add_vertices₂!(z_, 2)
add_edges!(z_, [1,2,2], [1,1,2])

line = UndirectedBipartiteGraph()
add_vertices₁!(line, 1)
add_vertices₂!(line, 2)
add_edges!(line, [1], [1])

parallel = UndirectedBipartiteGraph()
add_vertices₁!(parallel, 2)
add_vertices₂!(parallel, 2)
add_edges!(parallel, [1,2], [1,2])

merge = UndirectedBipartiteGraph()
add_vertices₁!(merge, 2)
add_vertices₂!(merge, 2)
add_edges!(merge, [1,2], [1,1])

Lspan = UndirectedBipartiteGraph()
add_vertices₁!(Lspan, 1)
add_vertices₂!(Lspan, 2)
add_edges!(Lspan, [1,1],[1,2])

I = UndirectedBipartiteGraph()
add_vertices₁!(I, 1)
add_vertices₂!(I, 2)

L = CSetTransformation(I, Lspan, V₁=[1], V₂=[1,2])
R = CSetTransformation(I, line, V₁=[1], V₂=[1,2])
m1 = CSetTransformation(Lspan, z_, V₁=[1], V₂=[1,2], E=[1, 2])
m2 = CSetTransformation(Lspan, z_, V₁=[1], V₂=[2,1], E=[2, 1])

@test is_isomorphic(parallel, rewrite_match(L, R, m1))
@test is_isomorphic(merge, rewrite_match(L, R, m2))

# Hypergraph homomorphisms
##########################

"""Theory of Acset schemas"""
@present TheoryACSetSchema(FreeSchema) begin
  (Name, AttrType) :: Data
  (HomOb, ObOb, AttrOb, DataOb, AIndex, HIndex) :: Ob
  (hsrc, htgt) :: Hom(HomOb,ObOb)
  asrc  :: Hom(AttrOb, ObOb)
  atgt  :: Hom(AttrOb,DataOb)
  hind  :: Hom(HIndex, HomOb)
  aind  :: Hom(AIndex, AttrOb)
  hname :: Attr(HomOb,Name)
  oname :: Attr(ObOb,Name)
  dname :: Attr(DataOb,Name)
  dtype :: Attr(DataOb,AttrType)
  aname :: Attr(AttrOb,Name)
end

const MetaCatDesc = ACSetType(TheoryACSetSchema, index=[
  :hsrc, :htgt, :asrc, :atgt, :hname, :oname, :dname, :aname]){Symbol, DataType}

"""Theory of Acset instances"""
@present TheoryACSet <: TheoryACSetSchema begin
  AttrVal :: Data
  (Row, FK, Datum) :: Ob
  rowob :: Hom(Row, ObOb)
  (fksrc, fktgt) :: Hom(FK, Row)
  fkhom :: Hom(FK, HomOb)
  drow :: Hom(Datum, Row)
  dattr :: Hom(Datum, AttrOb)
  dval :: Attr(Datum, AttrVal)
end

const MetaACSet = ACSetType(TheoryACSet, index=[
  :hsrc, :htgt, :asrc, :atgt, :hname, :oname, :dname, :aname]
  ){Symbol, DataType, Any}

"""Theory of Acset transformations"""
@present TheoryACSetTrans <: TheoryACSetSchema begin
  AttrVal :: Data

  # Src Data
  (Row1, FK1, Datum1) :: Ob
  rowob1 :: Hom(Row1, ObOb)
  (fksrc1, fktgt1) :: Hom(FK1, Row1)
  fkhom1 :: Hom(FK1, HomOb)
  drow1 :: Hom(Datum1, Row1)
  dattr1 :: Hom(Datum1, AttrOb)
  dval1 :: Attr(Datum1, AttrVal)

  # Tgt data
  (Row2, FK2, Datum2) :: Ob
  rowob2 :: Hom(Row2, ObOb)
  (fksrc2, fktgt2) :: Hom(FK2, Row2)
  fkhom2 :: Hom(FK2, HomOb)
  drow2 :: Hom(Datum2, Row2)
  dattr2 :: Hom(Datum2, AttrOb)
  dval2 :: Attr(Datum2, AttrVal)

  # Components
  comp :: Hom(Row1, Row2)
end

const MetaACSetTransformation = ACSetType(
  TheoryACSetTrans){Symbol, DataType, Any}

"""Convert an Acset description into its Acset representation"""
function to_cset(ad::AttrDesc, datatypes::Vector{Any}, inds::Vector{Symbol})::ACSet
  res = MetaCatDesc()
  catdesc, dtypes, anames, asrc, atgt = typeof(ad).parameters
  obj, homs, src, tgt = catdesc.parameters
  add_parts!(res, :ObOb, length(obj), oname=obj)
  add_parts!(res, :HomOb, length(homs), hname=homs, hsrc=src, htgt=tgt)
  add_parts!(res, :DataOb, length(dtypes), dname=dtypes,dtype=datatypes);
  add_parts!(res, :AttrOb, length(anames), aname=anames, asrc=asrc, atgt=atgt);
  hinds = filter(x-> x in homs, inds)
  ainds = filter(x -> x in anames, inds)
  add_parts!(res, :HIndex, length(hinds), hind=[findfirst(x->x==i, homs) for i in hinds]);
  add_parts!(res, :AIndex, length(ainds), aind=[findfirst(x->x==i, anames) for i in ainds]);
  return res
end

"""
For each Ob in the cset, copy its table into res
(only if the corresponding table in res is empty)
Optionally use a dict to map names from the first schema to the second

Helper function to construct MetaACSetTransformations
"""
function inject_cset(cset::ACSet, res::ACSet, mapping::Dict{Symbol, Symbol}=Dict{Symbol,Symbol}())::ACSet
  for (k, v) in zip(keys(cset.tables), cset.tables)
    k_ = get(mapping, k, k)
    if haskey(res.tables, k_) && isempty(res.tables[k_])
      add_parts!(res, k_, length(v))
    end
  end
  srchom, reshom = [typeof(x).parameters[2].parameters[1].parameters[2] for x in [cset, res]]
  for h in srchom
    h_ = get(mapping, h, h)
    if h_ in reshom
      set_subpart!(res, h_, cset[h])
    end
  end
  srcattr, resattr = [typeof(x).parameters[2].parameters[3] for x in [cset, res]]
  for a in srcattr
    a_ = get(mapping, a, a)
    if a_ in resattr
      set_subpart!(res, a_, cset[a])
    end
  end

  return res
end

"""
Create a CSet to describe the data of a CSet
"""
function to_cset(cset::ACSet)::Pair{ACSet, Dict{Pair{Symbol, Int}, Int}}
  tparams = typeof(cset).parameters;
  catdesc = to_cset(tparams[2](),
                    collect(tparams[3].parameters),
                    Vector{Symbol}(collect(tparams[4])));

  # Copy over info from the catdesc
  res = inject_cset(catdesc, MetaACSet())

  # map table+row pair to unique row index
  rowdict = Dict{Pair{Symbol, Int}, Int}()
  row_counter = 0
  for (k, v) in zip(keys(cset.tables), cset.tables)
    ob = catdesc.indices[:oname][k][1]
    for i in 1:length(v)
      row_counter+=1
      rowdict[k => i] = row_counter
      add_part!(res, :Row, rowob=ob)
    end
  end
  for (k, v) in zip(keys(cset.tables), cset.tables)
    ob = catdesc.indices[:oname][k][1]
    homs = catdesc.indices[:hsrc][ob]
    attrs = catdesc.indices[:asrc][ob]
    for (i, row) in enumerate(v)
      row_id = rowdict[k => i]
      row_counter += 1
      for morph in homs
        morphname = catdesc[:hname][morph]
        tgt_name  = catdesc[:oname][catdesc[:htgt][morph]]
        add_part!(res, :FK, fksrc=row_id, fkhom=morph,
                            fktgt=rowdict[tgt_name => row[morphname]])
      end
      for at in attrs
        attrname = catdesc[:aname][at]
        add_part!(res, :Datum, drow=row_id, dattr=at, dval=row[attrname])
      end
    end
  end

  return res => rowdict
end

"""
Create an ACSet to describe the data of a ACSetTransformation

Helper function to construct HΣrewrite
"""
function to_cset(m::ACSetTransformation)::Tuple{ACSet, Dict{Pair{Symbol, Int}, Int}, Dict{Pair{Symbol, Int}, Int}}
  tparams = typeof(dom(m)).parameters;
  catdesc = to_cset(tparams[2](),
                    collect(tparams[3].parameters),
                    Vector{Symbol}(collect(tparams[4])));

  # Copy over info from the catdesc
  res = inject_cset(catdesc, ACSetType(TheoryACSetTrans){Symbol, DataType, Any}())

  src, srcdict = to_cset(dom(m));
  tgt, tgtdict = to_cset(codom(m));

  allnames = ["Row", "FK", "Datum", "rowob", "fksrc", "fktgt", "fkhom", "drow", "dattr", "dval"]
  d1, d2 = [Dict([Symbol(x)=> Symbol(string(x)*i) for x in allnames]) for i in ["1", "2"]]
  res = inject_cset(src, res, d1)
  res = inject_cset(tgt, res, d2)

  # add components
  comps = repeat([0], nparts(res,:Row1))
  for (k, v) in zip(keys(components(m)), components(m))
    for (i, val) in enumerate(v.func)
      comps[srcdict[k=>i]] = tgtdict[k => val]
    end
  end
  set_subpart!(res, :comp, comps)
  return (res, srcdict, tgtdict)
end

"""
Given morphisms C1->C3<-C2, induce a morphism between the morphisms
via a morphism C1 -> C2

Helper function to construct HΣrewrite
"""
function transformation_of_transformations(m1::ACSetTransformation, m2::ACSetTransformation, msrc::ACSetTransformation)::ACSetTransformation
  (m1_, src1, tgt1), (m2_, src2, tgt2) = to_cset(m1), to_cset(m2)
  comps = Dict{Symbol, Vector{Int}}([comp => collect(1:nparts(m1_,comp))
                for comp in [:HomOb,:ObOb,:AttrOb, :DataOb, :AIndex, :HIndex]])
  r1, fk1, d1 = [repeat([0], nparts(m1_,x)) for x in [:Row1, :FK1, :Datum1]]
  r2, fk2, d2 = [collect(1:nparts(m2_,x)) for x in [:Row2, :FK2, :Datum2]]
  for (k, vs) in zip(keys(components(msrc)), components(msrc))
    for (i, v) in enumerate(vs.func)
      r1[src1[k => i]] = src2[k => v]
      # this will fail if msrc.domain has any FKs / attributes
      # for our example below, this is ok
    end
  end
  for (k, v) in zip([r1, fk1, d1, r2, fk2, d2],
                    [:Row1, :FK1, :Datum1, :Row2, :FK2, :Datum2])
    comps[v] = k
  end
  res = ACSetTransformation(m1_, m2_; comps...)
  return res
end

@present ThHypP(FreeSchema) begin
  (V, E₁₁, E₂₁)::Ob
  Name::Data
  (s, t)::Hom(E₁₁, V)
  (s₁, s₂, t₁)::Hom(E₂₁, V)
  Label₁::Attr(E₁₁, Name)
  Label₂::Attr(E₂₁, Name)
end

HypergraphP = ACSetType(ThHypP)
HΣ = @acset HypergraphP{Symbol} begin
  V = 1
  E₁₁ = 2
  E₂₁ = 1
  s = [1,1]
  t = [1,1]
  s₁ = [1]
  s₂ = [1]
  t₁ = [1]
  Label₁ = [:f, :g]
  Label₂ = [:plus]
end

H = @acset HypergraphP{Symbol} begin
  V = 9
  E₁₁ = 4
  E₂₁ = 2
  s = [1,2,3,6]
  t = [4,5,6,8]
  s₁ = [4,7]
  s₂ = [5,8]
  t₁ = [7,9]
  Label₁ = [:f, :g, :f, :g]
  Label₂ = [:plus,:plus]
end

ϕ = ACSetTransformation(H, HΣ, V=ones(Int, 9), E₁₁ = [1,2,1,2], E₂₁ = ones(Int,2))
(m1_,Hdict, HΣdict) = to_cset(ϕ);

HL = @acset HypergraphP{Symbol} begin
  V = 5
  E₂₁ = 2
  s₁ = [1,3]
  s₂ = [2,4]
  t₁ = [3,5]
  Label₂ = [:plus,:plus]
end

HK = @acset HypergraphP{Symbol} begin
  V = 3
end

HR = @acset HypergraphP{Symbol} begin
  V = 5
  E₂₁ = 2
  s₁ = [2,1]
  s₂ = [3,4]
  t₁ = [4,5]
  Label₂ = [:plus,:plus]
end

L = ACSetTransformation(HK, HL, V=[1,2,4])
R = ACSetTransformation(HK, HR, V=[1,2,3])
m= ACSetTransformation(HL, H, V=[4,5,7,8,9], E₂₁ = [1,2])
Hrewrite = rewrite_match(L,R,m); # WORKS
ϕHrewrite = ACSetTransformation(Hrewrite, HΣ, V=ones(Int,9), E₁₁ = [1,2,1,2], E₂₁ = [1,1])
ϕHrewrite_ = to_cset(ϕHrewrite)[1]; # convert to acset to check isomorphism later

ϕL = ACSetTransformation(HL, HΣ, V=ones(Int, 5),E₂₁ = [1,1]);
ϕK = ACSetTransformation(HK, HΣ, V=ones(Int, 3));
ϕR = ACSetTransformation(HR, HΣ, V=ones(Int, 5), E₂₁ = [1,1]);

(ϕ_,Hdict, HΣdict) = to_cset(ϕ);
(L_,HLdict, HLΣdict) = to_cset(ϕL);
(K_,HKdict, HKΣdict) = to_cset(ϕK);
(R_,HRdict, HRΣdict) = to_cset(ϕR);
l_ = transformation_of_transformations(ϕK, ϕL, L);
r_ = transformation_of_transformations(ϕK, ϕR, R);
m_ = homomorphism(L_, ϕ_); # there is only 1
HΣrewrite = rewrite_match(l_,r_,m_);

@test is_isomorphic(ϕHrewrite_, HΣrewrite)

end


text:
@test get_discretization_counts(LinearDiscretizer([0.0,1.0,2.0]), [0.5,0.5,0.5,1.5,1.5]) == [3,2]

# Test for multidimensional discretization
disc1 = LinearDiscretizer([0.0, 1.0, 2.0])
disc2 = LinearDiscretizer([2.0, 3.0, 4.0, 5.0])
data = [0.5 2.5
        0.5 2.5
        0.5 3.5
        0.5 4.5
        1.5 2.5
        1.5 3.5
        1.5 3.5
        1.5 3.5]
counts = [2 1 1
          1 3 0]
@test get_discretization_counts([disc1, disc2], data) == counts


text:
"""
Optical properties abstract struct
"""
abstract type OpticalProperties <: AbstractInterceptionModel end


text:
export z, t, abbott

include("end_effects.jl")

"""
    Rz(R, λc=250)

Determine roughness of profile-component `component`. Filterwavelength in accordance with [DIN EN ISO 3274:1998](DIN_EN_ISO_3274_1998.pdf).
"""
function z(profile::SurfaceProfile, component::Symbol, eliminate_end_effect=zeropad)
    c = eliminate_end_effect(profile, component)
    bc = profile.λc * profile.resolution
    maximum(-[-(extrema(c[n+1:n+bc])...) for n ∈ 1:bc:(length(c)-bc)])
end

"""
    Wt(W, λc=250)

Determine waviness of the waviness-profile W. Filterwavelength in accordance with [DIN EN ISO 3274:1998](DIN_EN_ISO_3274_1998.pdf).
"""
function t(profile::SurfaceProfile, component::Symbol, eliminate_end_effect=zeropad)
    c = eliminate_end_effect(profile, component)
    maximum(c) - minimum(c)
end

"""
    abbott(surf)

Calculate an [Abbott-Firestone-Curve](https://en.wikipedia.org/wiki/Abbott-Firestone_curve) for the line scan given in `surf`.
"""
abbott(surf) = [quantile(surf, i) for i ∈ 1:-0.01:0]


text:
export RandomizedKaczmarz, RandomizedBlockKaczmarz

"""
Randomized Kaczmarz

Section 3.3 in Gower, R. M., & Richtárik, P. (2015).
Randomized iterative methods for linear systems.
SIAM Journal on Matrix Analysis and Applications, 36(4), 1660-1690.

RK takes a step in the direction of the negative stochastic gradient.
This means that it is equivalent to the SGD method.
However, the stepsize choice is very special: RK chooses the stepsize which leads to the point
which is closest to x* in the Euclidean norm.
"""
function RandomizedKaczmarz(stp::AbstractStopping; is_zero_start::Bool = false, kwargs...)
  A, b = get_matrix(stp.pb), get_vector(stp.pb)
  state = stp.current_state
  m = size(A, 1)
  if is_zero_start
    state.res .= b
  else
    state.res .= b .- A * state.x
  end
  OK = start!(stp, no_opt_check = true)

  while !OK
    i = Int(floor(rand() * m) + 1)
    Ai = view(A, i, :)
    AiAi = dot(Ai, Ai)
    if AiAi != 0
      Aix = dot(Ai, state.x)
      state.x .-= (Aix - b[i]) / AiAi * Ai
    end
    state.res = b - A * state.x
    OK = cheap_stop!(stp)
  end

  return stp
end

"""
Randomized block Kaczmarz

Section 3.5 in Gower, R. M., & Richtárik, P. (2015).
Randomized iterative methods for linear systems.
SIAM Journal on Matrix Analysis and Applications, 36(4), 1660-1690.
"""
function RandomizedBlockKaczmarz(
  stp::AbstractStopping;
  is_zero_start::Bool = false,
  r::Int = 15,
  rand_r::Bool = false,
  kwargs...,
)
  A, b = get_matrix(stp.pb), get_vector(stp.pb)
  state = stp.current_state
  m = size(A, 1)
  if is_zero_start
    state.res .= b
  else
    state.res .= b .- A * state.x
  end
  OK = start!(stp, no_opt_check = true)

  while !OK
    r = rand_r ? Int(floor(rand() * m) + 1) : min(r, m)
    sub = zeros(Int64, r)
    #sample!(1:m, sub) #x is a subset of [1,...,m] of size r #sample repeat entries
    sub = StatsBase.randperm(m)[1:r]

    Ai = view(A, sub, :)
    resk = state.res[sub]
    if !(Ai == 0)
      state.x .-= Ai' * pinv(Matrix(Ai * Ai')) * resk
    end
    state.res .= A * state.x .- b
    OK = cheap_stop!(stp)
  end

  return stp
end


text:
struct FEOperatorFromWeakForm <: FEOperator
  res::Function
  jac::Function
  trial::FESpace
  test::FESpace
  assem::Assembler
  function FEOperatorFromWeakForm(
    res::Function,
    jac::Function,
    trial::FESpace,
    test::FESpace,
    assem::Assembler)
    @assert ! isa(test,TrialFESpace) """\n
    It is not allowed to build a FEOperator with a test space of type TrialFESpace.

    Make sure that you are writing first the trial space and then the test space when
    building an AffineFEOperator or a FEOperatorFromWeakForm.
    """
    new(res,jac,trial,test,assem)
  end
end

function FEOperator(
  res::Function,jac::Function,trial::FESpace,test::FESpace,assem::Assembler)
  FEOperatorFromWeakForm(res,jac,trial,test,assem)
end

function FEOperator(res::Function,jac::Function,args...)
  assem = SparseMatrixAssembler(args...)
  trial,test, = args
  FEOperator(res,jac,trial,test,assem)
end

function FEOperator(
  res::Function,trial::FESpace,test::FESpace,assem::Assembler)
  jac(u,du,dv) = jacobian(x->res(x,dv),u)
  FEOperatorFromWeakForm(res,jac,trial,test,assem)
end

function FEOperator(res::Function,args...)
  assem = SparseMatrixAssembler(args...)
  trial,test, = args
  FEOperator(res,trial,test,assem)
end

get_test(op::FEOperatorFromWeakForm) = op.test
get_trial(op::FEOperatorFromWeakForm) = op.trial

function allocate_residual(op::FEOperatorFromWeakForm,uh::FEFunction)
  V = get_test(op)
  v = get_fe_basis(V)
  vecdata = collect_cell_vector(V,op.res(uh,v))
  allocate_vector(op.assem, vecdata)
end

function residual!(b::AbstractVector,op::FEOperatorFromWeakForm,uh::FEFunction)
  V = get_test(op)
  v = get_fe_basis(V)
  vecdata = collect_cell_vector(V,op.res(uh,v))
  assemble_vector!(b,op.assem, vecdata)
  b
end

function allocate_jacobian(op::FEOperatorFromWeakForm,uh::FEFunction)
  U = get_trial(op)
  V = get_test(op)
  du = get_trial_fe_basis(U)
  v = get_fe_basis(V)
  matdata = collect_cell_matrix(U,V,op.jac(uh,du,v))
  allocate_matrix(op.assem, matdata)
end

function jacobian!(A::AbstractMatrix,op::FEOperatorFromWeakForm,uh::FEFunction)
  U = get_trial(op)
  V = get_test(op)
  du = get_trial_fe_basis(U)
  v = get_fe_basis(V)
  matdata = collect_cell_matrix(U,V,op.jac(uh,du,v))
  assemble_matrix!(A,op.assem,matdata)
  A
end

function residual_and_jacobian!(
  b::AbstractVector,A::AbstractMatrix,op::FEOperatorFromWeakForm,uh::FEFunction)
  U = get_trial(op)
  V = get_test(op)
  du = get_trial_fe_basis(U)
  v = get_fe_basis(V)
  data = collect_cell_matrix_and_vector(U,V,op.jac(uh,du,v),op.res(uh,v))
  assemble_matrix_and_vector!(A, b, op.assem, data)
  (b,A)
end

function residual_and_jacobian(op::FEOperatorFromWeakForm,uh::FEFunction)
  U = get_trial(op)
  V = get_test(op)
  du = get_trial_fe_basis(U)
  v = get_fe_basis(V)
  data = collect_cell_matrix_and_vector(U,V,op.jac(uh,du,v),op.res(uh,v))
  A, b = assemble_matrix_and_vector(op.assem, data)
  (b, A)
end


text:
module BigramMean

require("Corpora.jl")
require("CommonTools.jl")
require("CVMAbstract.jl")
require("DataTypes.jl")

## Imports
using DataTypes
using Base.LinAlg.BLAS.axpy!
using Corpora.Token
using Corpora.text
using CommonTools: tanh, tanh!, ultra_fast_tanh, ultra_fast_tanh!
using CommonTools: hard_tanh, hard_tanh!
using CommonTools.CVMIntermediates

using CVMAbstractModel
import CVMAbstractModel: compose!, store_gradients!

## Exports
export BigramMeanModel
export BigramMeanParameters
export compose, step!, compose_w_l2
export compose!
export store_gradients!

immutable BigramMeanModel <: AbstractModel
    wordvecdim::Int64;
    margin::Float;
    regularizer::Float;
    lr_init::Float;
    lr_min::Float;
    niter::Int64;
    gamma::Float;
    nnoise::Int64;
    batchsize::Int64;
end

# Stores intermediate values
type BigramMeanParameters <: AbstractParameters
    wordvecs::Matrix{Float};
    gradients::Matrix{Float};
    rms::Matrix{Float};
    pivot_vecs_l::Matrix{Float};
    pivot_gradients_l::Matrix{Float};
    pivot_rms_l::Matrix{Float};
    pivot_vecs_r::Matrix{Float};
    pivot_gradients_r::Matrix{Float};
    pivot_rms_r::Matrix{Float};
end

function BigramMeanParameters(wordvecs::Matrix{Float},
                            gradients::Matrix{Float},
                            rms::Matrix{Float})
    return BigramMeanParameters(wordvecs, gradients, rms,
                              Array(Float,0,0), Array(Float,0,0),
                              Array(Float,0,0),
                              Array(Float,0,0), Array(Float,0,0),
                              Array(Float,0,0))
end

BigramMeanParameters() = BigramMeanParameters(Array(Float,0,0), Array(Float,0,0),
                                                         Array(Float,0,0))

function constrastive_update1!(o::Matrix{Float64}, idx::Int64,
                                 pre_grad::Vector{Float64},
                                 target::DenseArray{Float64,1},
                                 scaling::Float64)
    """
        Propagates the given ``pre_grad`` gradient through the first
        or last vector bigram of the a sentence.
    """
    @simd for i=1:size(target,1)
        # We can compute the derivative of of tanh(x) = o as 1-o^2
        # Since o is elementwise squared already, it reduces to (1-o)
        # o gets elementwise square in the function that calls this method!
        target[i] += (1.0 - o[i,idx]) * pre_grad[i] * scaling;
    end
end

function constrastive_update2!(o::Matrix{Float64}, idx1::Int64, idx2::Int64,
                                 pre_grad::Vector{Float64},
                                 target::DenseArray{Float64,1},
                                 scaling::Float64)
    """
        Propagates the given ``pre_grad`` gradient through a vector bigram
        in the middle of a sentence.
    """
    @simd for i=1:size(target,1)
        target[i] += (2.0 - o[i,idx1] - o[i,idx2]) * pre_grad[i] * scaling;
    end
end

function compose_w_l2(model::BigramMeanModel, s::Vector{Token},
                        o::Matrix{Float64}, meanvec::Vector{Float64})
    """
        Compose a vector representation for a sentence while calculating
        the l2 norm of the sentence's word vectors at the same time.
        The sentence is given as a list of tokens.
    """
    l2reg = l2(s[1].vector)
    for iword=2:size(s,1)
        iw1 = iword-1
        v1 = s[iword-1].vector
        v2 = s[iword].vector
        @simd for io=1:model.wordvecdim
            @inbounds o[io,iw1] = tanh(v1[io] + v2[io])
        end
        #l2reg += sum(v2.^2)
        l2reg += l2(v2)
    end

    # Get column mean of o
    size_o_2 = size(s,1)-1
    size_o_1 = model.wordvecdim
    Base.LinAlg.BLAS.scal!(model.wordvecdim,0.0,meanvec,1);
    for y=1:size_o_2
        @simd for x=1:size_o_1
            @inbounds meanvec[x] += o[x,y]
        end
    end
    return meanvec, o, l2reg
end

function compose_w_l2(model::BigramMeanModel, s::Vector{Token},
                        o::Matrix{Float64})
    return compose_w_l2(model, s, o, zeros(Float64, model.wordvecdim));
end

function compose_w_l2(model::BigramMeanModel, s::Vector{Token})
    return compose_w_l2(model, s, zeros(Float64, model.wordvecdim, size(s,1)-1));
end

function addToColTanh!(icol, v1, v2, o)
    """
        High-speed function that sums up two vectors and
        applies tanh componentwise to the result.
    """
    @simd for i=1:length(v1)
        @inbounds o[i, icol] = v1[i]
        @inbounds o[i, icol] += v2[i]
        @inbounds o[i, icol] = ultra_fast_tanh(o[i,icol])
    end
end

function compose!(model::BigramMeanModel, s::Vector{Token},
                                                    o::Matrix{Float64},
                                                    result::Vector{Float64})
    """
        Compose a sentence vector from a list of tokens while
        storing intermediates results.
    """
    # BOTH: o and result are modified in this function

    # Create the intermediary result matrix "o"
    for iword=2:size(s,1)
        iw1 = iword-1
        v1 = s[iw1].vector
        v2 = s[iword].vector
        addToColTanh!(iw1, v1, v2, o)
    end

    # Add the columns of o up to get the final result
    size_o_2 = size(s,1)-1
    size_o_1 = model.wordvecdim
    Base.LinAlg.BLAS.scal!(model.wordvecdim,0.0,result,1);
    for y=1:size_o_2
        @simd for x=1:size_o_1
            result[x] += o[x,y]
        end
    end
    Base.LinAlg.BLAS.scal!(model.wordvecdim,1.0/size_o_2,result,1);
    return result
end

function compose(model::BigramMeanModel, s::Vector{Token})
    return compose!(model, s, zeros(model.wordvecdim, size(s,1)),
                    zeros(Float64, model.wordvecdim))
end

function compose!(model::BigramMeanModel, s::Vector{Token},
                                    intermediate::CVMIntermediates)
    return compose!(model, s, intermediate.ofull, intermediate.result)
end

function store_gradients!(model::BigramMeanModel, pre_grad::Vector{Float64},
                    s::Vector{Token}, intermediate::CVMIntermediates;
                    intermediates_precomputed::Bool=false)
    """
        Given a gradient from a higher level in the neural network,
        back-propagate the gradient through to the word representations.
    """

    #o = copy(intermediate.ofull)
    o = intermediate.ofull
    size_o = size(s,1)-1
    # This line makes the difference between addition and mean
    scaling = 1.0/size_o

    if !intermediates_precomputed
        # Square the output (because tanh'(x) = 1 - tanh(x)^2)
        for y=1:size_o
            @simd for x=1:size(o,1)
                @inbounds o[x,y] *= o[x,y]
            end
        end
    end

    constrastive_update1!(o, 1, pre_grad, s[1].gradient,
                            scaling);
    for iword=2:size_o
        constrastive_update2!(o, iword-1, iword, pre_grad, s[iword].gradient,
                                scaling)
    end
    constrastive_update1!(o, size_o, pre_grad, s[end].gradient,
                            scaling);
end

end

max_stars_repo_name: ogh/binclusion

max_issues_repo_name: ogh/binclusion

max_forks_repo_name: ogh/binclusion


text:
using GFF3
using Test
using Documenter

using BioGenerics
using FASTX.FASTA
using FormatSpecimens
using GenomicFeatures

import BioSequences: @dna_str

import BGZFStreams

import BioGenerics.Exceptions: MissingFieldException

@testset "GFF3" begin
    record = GFF3.Record()
    @test !BioGenerics.isfilled(record)
    @test repr(record) == "GFF3.Record: <not filled>"

    record = GFF3.Record("CCDS1.1\tCCDS\tgene\t801943\t802434\t.\t-\t.\tNAME=LINC00115")
    @test BioGenerics.isfilled(record)
    @test GFF3.isfeature(record)
    @test BioGenerics.hasseqname(record)
    @test GFF3.hasseqid(record)
    @test BioGenerics.seqname(record) == GFF3.seqid(record) == "CCDS1.1"
    @test GFF3.hassource(record)
    @test GFF3.source(record) == "CCDS"
    @test GFF3.hasfeaturetype(record)
    @test GFF3.featuretype(record) == "gene"
    @test GFF3.hasseqstart(record) === BioGenerics.hasleftposition(record) === true
    @test GFF3.seqstart(record) === BioGenerics.leftposition(record) === 801943
    @test GFF3.hasseqend(record) === BioGenerics.hasrightposition(record) === true
    @test GFF3.seqend(record) === BioGenerics.rightposition(record) === 802434
    @test !GFF3.hasscore(record)
    @test_throws MissingFieldException GFF3.score(record)
    @test GFF3.hasstrand(record)
    @test strand(record) === GFF3.strand(record) === STRAND_NEG
    @test !GFF3.hasphase(record)
    @test_throws MissingFieldException GFF3.phase(record)
    @test GFF3.attributes(record) == ["NAME" => ["LINC00115"]]
    @test GFF3.attributes(record, "NAME") == ["LINC00115"]
    @test GFF3.content(record) == "CCDS1.1\tCCDS\tgene\t801943\t802434\t.\t-\t.\tNAME=LINC00115"
    @test startswith(repr(record), "GFF3.Record:\n")
    @test string(record) == "CCDS1.1\tCCDS\tgene\t801943\t802434\t.\t-\t.\tNAME=LINC00115"

    record = GFF3.Record("##gff-version 3")
    @test BioGenerics.isfilled(record)
    @test GFF3.isdirective(record)
    @test GFF3.content(record) == "gff-version 3"
    @test convert(String, record) == "##gff-version 3"

    record = GFF3.Record("#comment")
    @test BioGenerics.isfilled(record)
    @test GFF3.iscomment(record)
    @test GFF3.content(record) == "comment"
    @test convert(String, record) == "#comment"

    function check_gff3_parse(filename)
        # Reading from a stream
        num_intervals = 0
        for interval in GFF3.Reader(open(filename))
            num_intervals += 1
        end

        # Reading from a regular file
        for interval in open(GFF3.Reader, filename)
        end

        collection = IntervalCollection(open(GFF3.Reader, filename))
        @test length(collection) == num_intervals

        # in-place parsing
        stream = open(GFF3.Reader, filename)
        entry = eltype(stream)()
        while !eof(stream)
            try
                read!(stream, entry)
            catch ex
                if isa(ex, EOFError)
                    break
                end
            end
        end
        close(stream)

        # copy
        records = GFF3.Record[]
        reader = open(GFF3.Reader, filename)
        output = IOBuffer()
        writer = GFF3.Writer(output)
        for record in reader
            write(writer, record)
            push!(records, record)
        end
        close(reader)
        flush(writer)

        records2 = GFF3.Record[]
        for record in GFF3.Reader(IOBuffer(take!(output)))
            push!(records2, record)
        end
        return records == records2
    end

    dir_gff3 = path_of_format("GFF3")

    for specimen in list_valid_specimens("GFF3")

        if hastag(specimen, "gzip")
            # skip compressed files
            continue
        end

        filepath = joinpath(dir_gff3, filename(specimen))

        @test check_gff3_parse(filepath)

    end

    for specimen in list_invalid_specimens("GFF3")

        if "gzip" ∈ split(get(specimen, "tags", ""))
            # skip compressed files
            continue
        end

        filepath = joinpath(dir_gff3, filename(specimen))

        @test_throws Exception check_gff3_parse(filepath)
    end

    # no fasta
    test_input = """
1	havana	exon	870086	870201	.	-	.	Parent=transcript:ENST00000432963;Name=ENSE00001791782;constitutive=0;ensembl_end_phase=-1;ensembl_phase=-1;exon_id=ENSE00001791782;rank=1;version=1
1	havana	lincRNA	868403	876802	.	-	.	ID=transcript:ENST00000427857;Parent=gene:ENSG00000230368;Name=FAM41C-002;biotype=lincRNA;havana_transcript=OTTHUMT00000007022;havana_version=1;transcript_id=ENST00000427857;transcript_support_level=3;version=1
"""
    stream = GFF3.Reader(IOBuffer(test_input))
    collect(stream)
    @test !GFF3.hasfasta(stream)
    @test_throws Exception GFF3.getfasta(stream)

    # implicit fasta
    test_input2 = string(test_input, """
>seq1
ACGTACGT
>seq2
TGCATGCA
""")
    stream = GFF3.Reader(IOBuffer(test_input2))
    collect(stream)
    @test GFF3.hasfasta(stream)
    @test collect(GFF3.getfasta(stream)) ==
        [FASTA.Record("seq1", dna"ACGTACGT")
         FASTA.Record("seq2", dna"TGCATGCA")]

    # explicit fasta
    test_input3 = string(test_input, """
##FASTA
>seq1
ACGTACGT
>seq2
TGCATGCA
""")
    stream = GFF3.Reader(IOBuffer(test_input3))
    collect(stream)
    @test GFF3.hasfasta(stream)
    @test collect(GFF3.getfasta(stream)) ==
        [FASTA.Record("seq1", dna"ACGTACGT")
         FASTA.Record("seq2", dna"TGCATGCA")]

    test_input4 = """
##directive1
#comment1
##directive2
1	havana	exon	869528	869575	.	-	.	Parent=transcript:ENST00000432963;Name=ENSE00001605362;constitutive=0;ensembl_end_phase=-1;ensembl_phase=-1;exon_id=ENSE00001605362;rank=2;version=1
1	havana	exon	870086	870201	.	-	.	Parent=transcript:ENST00000432963;Name=ENSE00001791782;constitutive=0;ensembl_end_phase=-1;ensembl_phase=-1;exon_id=ENSE00001791782;rank=1;version=1
##directive3
#comment2
##directive4
1	havana	lincRNA	868403	876802	.	-	.	ID=transcript:ENST00000427857;Parent=gene:ENSG00000230368;Name=FAM41C-002;biotype=lincRNA;havana_transcript=OTTHUMT00000007022;havana_version=1;transcript_id=ENST00000427857;transcript_support_level=3;version=1
##directive5
#comment3
##directive6
"""
    stream = GFF3.Reader(IOBuffer(test_input4), save_directives=true)
    read(stream)
    @test GFF3.directives(stream) == ["directive1", "directive2"]
    read(stream)
    @test isempty(GFF3.directives(stream))
    read(stream)
    @test GFF3.directives(stream) == ["directive3", "directive4"]
    @test_throws EOFError read(stream)
    @test eof(stream)
    close(stream)
    @test GFF3.directives(stream) == ["directive5", "directive6"]

    test_input5 = """
    ##directive1
    feature1\t.\t.\t.\t.\t.\t.\t.\t
    #comment1
    feature2\t.\t.\t.\t.\t.\t.\t.\t
    ##directive2
    feature3\t.\t.\t.\t.\t.\t.\t.\t
    """
    @test [r.kind for r in GFF3.Reader(IOBuffer(test_input5))] == [:feature, :feature, :feature]
    @test [r.kind for r in GFF3.Reader(IOBuffer(test_input5), skip_directives=false)] == [:directive, :feature, :feature, :directive, :feature]
    @test [r.kind for r in GFF3.Reader(IOBuffer(test_input5), skip_directives=false, skip_comments=false)] == [:directive, :feature, :comment, :feature, :directive, :feature]

    @testset "eachoverlap" begin
        path = joinpath(path_of_format("GFF3"), "TAIR10.part.gff.bgz")
        stream = BGZFStreams.BGZFStream(path)
        reader = GFF3.Reader(stream, index=string(path, ".tbi"))
        for (interval, n_records) in [
                (Interval("ChrC", 1:70_000), 68),
                (Interval("Chr4", 1:10_000),  2),
                (Interval("ChrM", 1:30_000), 10),]
            n = 0
            for record in eachoverlap(reader, interval)
                n += 1
            end
            @test n == n_records
        end
        @test isa(GFF3.Reader(path), GFF3.Reader)
    end

    # Include doctests.
    DocMeta.setdocmeta!(GFF3, :DocTestSetup, :(using GFF3); recursive=true)
    doctest(GFF3; manual = false)
end


text:
"""
# Method p46

Stability (buckling) analysis of elastic beams using a 2-node Beam
structural element and Line finite elements. Elastic foundation is
optional.

### Constructors
```julia
p46(data)
```
### Arguments
```julia
* `data::Dict{Symbol, Any}` : Dictionary containing all input data
```

### Required data dictionary keys
```julia
* struc_el::StructuralElement                          : Type of  structural fin_el
* support::Array{Tuple{Int,Array{Int,1}},1}        : Fixed-displacements vector
* properties::Vector{Float64}                          : Material properties
* x_coords::FloatRange{Float64}                        : x-coordinate vector
```

### Optional additional data dictionary keys
```julia
* etype::Vector{Int}         : Element material vector if np_types > 1
* limit = 250                  : Iteration limit
* tol = 0.0001                 : Tolerance for iteration convergence
```

### Related help
```julia
?StructuralElement             : List of available structural element types
?Beam                          : Help on a Beam structural element
?FiniteElement                 : List finite element types
?Line                          : Help on Line finite element
```
"""
function p46(data::Dict{Symbol, Any})

  # Parse & check FEdict data

  if :struc_el in keys(data)
    struc_el = data[:struc_el]
  else
    println("No fin_el type specified.")
    return
  end

  ndim = struc_el.ndim
  nst = struc_el.nst

  # Add radial stress
  if ndim == 3 && struc_el.axisymmetric
    nst = 4
  end

  fin_el = struc_el.fin_el
  @assert typeof(fin_el) <: FiniteElement

  if typeof(fin_el) == Line
    (nels, nn) = mesh_size(fin_el, struc_el.nxe)
  elseif typeof(fin_el) == Triangle || typeof(fin_el) == Quadrilateral
    (nels, nn) = mesh_size(fin_el, struc_el.nxe, struc_el.nye)
  elseif typeof(fin_el) == Hexahedron
    (nels, nn) = mesh_size(fin_el, struc_el.nxe, struc_el.nye, struc_el.nze)
  else
    println("$(typeof(fin_el)) is not a known finite element.")
    return
  end

  nodof = ndim                    # Degrees of freedom per node
  ndof = fin_el.nod * nodof      # Degrees of freedom per fin_el

  # Update penalty if specified in FEdict

  penalty = 1e20
  if :penalty in keys(data)
    penalty = data[:penalty]
  end

  # Allocate all arrays

  # Start with arrays to be initialized from FEdict

  if :properties in keys(data)
    prop = zeros(size(data[:properties], 1), size(data[:properties], 2))
    for i in 1:size(data[:properties], 1)
      prop[i, :] = data[:properties][i, :]
    end
  else
    println("No :properties key found in FEdict")
  end

  nf = ones(Int, nodof, nn)
  if :support in keys(data)
    for i in 1:size(data[:support], 1)
      nf[:, data[:support][i][1]] = data[:support][i][2]
    end
  end

  x_coords = zeros(nn)
  if :x_coords in keys(data)
    x_coords = data[:x_coords]
  end

  y_coords = zeros(nn)
  if :y_coords in keys(data)
    y_coords = data[:y_coords]
  end

  z_coords = zeros(nn)
  if :z_coords in keys(data)
    z_coords = data[:z_coords]
  end

  etype = ones(Int, nels)
  if :etype in keys(data)
    etype = data[:etype]
  end

  g_num = zeros(Int, fin_el.nod, nels)
  if :g_num in keys(data)
    g_num = data[:g_num]
  end

  # All other arrays

  points = zeros(struc_el.nip, ndim)
  g = zeros(Int, ndof)
  g_coord = zeros(ndim,nn)
  fun = zeros(fin_el.nod)
  coord = zeros(fin_el.nod, ndim)
  gamma = zeros(nels)
  jac = zeros(ndim, ndim)
  der = zeros(ndim, fin_el.nod)
  deriv = zeros(ndim, fin_el.nod)
  bee = zeros(nst,ndof)
  km = zeros(ndof, ndof)
  mm = zeros(ndof, ndof)
  gm = zeros(ndof, ndof)
  kg = zeros(ndof, ndof)
  eld = zeros(ndof)
  weights = zeros(struc_el.nip)
  g_g = zeros(Int, ndof, nels)
  num = zeros(Int, fin_el.nod)
  actions = zeros(ndof, nels)
  displacements = zeros(size(nf, 1), ndim)
  gc = ones(ndim, ndim)
  dee = zeros(nst,nst)
  sigma = zeros(nst)
  axial = zeros(nels)

  # Set global coordinates

  formnf!(nodof, nn, nf)
  neq = maximum(nf)
  evec = zeros(neq+1)

  ell = zeros(nels)
  if :x_coords in keys(data)
    for i in 1:length(data[:x_coords])-1
      ell[i] = data[:x_coords][i+1] - data[:x_coords][i]
    end
  end

  for i in 1:nels
    num = [i; i+1]
    num_to_g!(num, nf, g)
    g_g[:, i] = g
  end

  println("There are $(neq) equations.\n")

  gsm = spzeros(neq, neq)
  ggm = spzeros(neq, neq)
  for i in 1:nels
    km = beam_km(prop[etype[i], 1], ell[i])
    if size(prop, 2) > 1
      mm = beam_mm(prop[etype[i], 2], ell[i])
    end
    kg = beam_gm(ell[i])
    g = g_g[:, i]
    fsparm!(gsm, g, km+mm)
    fsparm!(ggm, g, kg)
  end

  limit = 10
  if :limit in keys(data)
    limit = data[:limit]
  end

  tol = 1.0e-5
  if :tol in keys(data)
    tol = data[:tol]
  end

  # No updates to arguments, drop the ! mark.
  (iters, evec, ival) = stability(gsm, ggm, tol, limit)
  evec[1] = 0.0
  (ival, iters, nn, evec, nf)
end

max_stars_repo_path: src/4 Static Equilibrium/p46.jl

max_issues_repo_path: src/4 Static Equilibrium/p46.jl

max_forks_repo_path: src/4 Static Equilibrium/p46.jl


text:
using JLD
using SharedArrays
using Distributed
include("../examples/baker.jl")
function obj_fun(x,y)
		return cos(4*y)
end
function obj_fun_erg_avg(s)
	nSteps = 10000
	u = 2*pi*rand(2)
	J = 0.
	u_trj = step(u,s,nSteps-1)
	x, y = view(u_trj,1,:), view(u_trj,2,:)
	J = sum(obj_fun.(x,y)/nSteps)
	return J
end
function get_Javg_vs_s(ind)
	s = zeros(4)
	n_pts = 10
	n_rep = 16000
	s_ind = LinRange(0.,1.0,n_pts)
	J = zeros(n_pts)
	J_proc = SharedArray{Float64}(n_rep)
	for i = 1:n_pts
		@show s_ind[i]
		s[ind] = s_ind[i]
		#s[3] = s_ind[i]
		J_proc .= 0.
		t = @distributed for n=1:n_rep
			J_proc[n] = obj_fun_erg_avg(s)/n_rep
		end
		wait(t)
		J[i] = sum(J_proc)
	end
	save("../data/obj_erg_avg/cos4y_s$(ind)_sens.jld",
		 "s$ind", s_ind,
		"J", J)
end

max_stars_repo_name: nishaChandramoorthy/linearResponse

max_issues_repo_name: nishaChandramoorthy/linearResponse

max_forks_repo_name: nishaChandramoorthy/linearResponse