text:
#!/usr/bin/env bash
    #=
    exec julia --project="$(realpath $(dirname $(realpath $(dirname $0))))" --color=yes --startup-file=no -e "include(popfirst!(ARGS))" \
    "${BASH_SOURCE[0]}" "$@"
    =#

include(joinpath(dirname(dirname(@__FILE__)), "src", "CodingTheory.jl"))

using .CodingTheory
import .CodingTheory.deepsym, .CodingTheory.displaymatrix

using Test
using Polynomials

function aresamearrays(A::AbstractVector{T}, B::AbstractVector{R}) where {T, R}
	A === B && return true
	length(A) ≠ length(B) && return false

	for (a, b) in zip(sort(A), sort(B))
		if a ≠ b
			return false
		end
	end

	return true
end

"""
Given two arrays A and B, ensures they have the same elements
"""
function aresamearrays(A::AbstractArray{T, N}, B::AbstractArray{R, M}) where {T, R, N, M}
	A === B && return true
	length(A) ≠ length(B) && return false

	for (a, b) in zip(sort(reshape(A, :)), sort(reshape(B, :)))
	# for (a, b) in zip(sort(A, dims = 1), sort(B, dims = 1))
		if a ≠ b
			return false
		end
	end

	return true
end

@time @testset "CodingTheory.jl" begin
	@test has_identity([1 0 0 2 3 0; 0 1 0 1 2 2; 0 0 1 4 3 0], 1) == true
	@test has_identity([1 0 0 2 3 0; 0 1 0 1 2 2; 0 0 1 4 3 0]) == true
	@test has_identity([1 0 1 2 3 0; 0 1 0 1 2 2; 0 0 1 4 3 0]) == true
	@test has_identity([1 0 1 2 3 0; 0 1 0 1 2 2; 0 0 1 4 3 0], 3) == false
	@test has_identity([-96 -66 20 1 0 0; -65 59 -82 0 1 0; -16 87 -113 0 0 1]) == true
	@test has_identity([78 -99 125 -123 -111 -71 17; -115 78 40 -88 81 -40 78; -99 126 -54 1 0 0 24; -55 88 42 0 1 0 -8; 119 55 2 0 0 1 -92; -40 -21 -89 -79 59 -44 9]) == true
	@test has_identity([78 -99 125 -123 -111 -71 17; -115 78 40 -88 81 -40 78; -99 126 -54 1 0 0 24; -55 88 42 0 1 0 -8; 119 55 2 0 0 1 -92; -40 -21 -89 -79 59 -44 9], 4) == false
	@test has_identity_on_left([1 0 0 2 3 0; 0 1 0 1 2 2; 0 0 1 4 3 0]) == true
	@test has_identity_on_left([-96 -66 20 1 0 0; -65 59 -82 0 1 0; -16 87 -113 0 0 1]) == false

	@test hamming_distance("ABC", "DBC") == 1
	@test hamming_distance("ABC", "DEF") == 3

	@test aresamearrays(hamming_ball([[1, 0, 1], [0, 1, 1], [1, 0, 0]], [1, 0, 0], 2), deepsym([[1, 0, 1], [1, 0, 0]]))
	@test aresamearrays(hamming_ball(list_span([2, 2, 2], [1, 2, 0], [0, 2, 1], 3), [1, 0, 0], 3), deepsym([[0, 0, 0], [0, 2, 1], [0, 1, 2], [1, 2, 0], [1, 1, 1], [1, 0, 2], [2, 1, 0], [2, 0, 1], [2, 2, 2]]))

	@test rate(3, 5, 4) ≈ 0.3662433802

	@test t_error_correcting([[0, 0, 0, 0], [2, 2, 2, 2]], 1) == true
	@test t_error_correcting([[1, 0, 1], [0, 1, 1], [1, 0, 0], [1, 1, 1]], 3) == false
	@test t_error_detecting([[1, 0, 1], [0, 1, 1], [1, 0, 0], [1, 1, 1]], 3) == false

	@test find_error_detection_max([[0, 0, 0, 0], [0, 1, 1, 1], [1, 0, 1, 0], [1, 1, 0, 1]], 2) == 1
	@test find_error_correction_max([[0, 0, 0, 0], [0, 1, 1, 1], [1, 0, 1, 0], [1, 1, 0, 1]], 2) == 0
	@test find_error_correction_max(list_span([1, 0, 1, 0], [0, 1, 1, 1], 2), 2) == 0
	@test find_error_detection_max(list_span([1, 0, 1, 0], [0, 1, 1, 1], 2), 2) == 1

	@test Polynomial([1, 2, 3, 4, 5, 6, 7, 8, 9], 3) == Polynomial([1, 2, 0, 1, 2, 0, 1, 2])

	@test isirreducible(Polynomial([1, 1, 0, 0, 1]), 2) == true
	@test isirreducible(Polynomial([1, 1, 1, 0, 1, 1]), 2) == true
	@test isirreducible(Polynomial([4, 4, 1]), 2) == false
	@test isirreducible(Polynomial([1, 0, 1]), 2) == false
	@test isirreducible(Polynomial([-2, 0, 1]), 2) == false
	@test isirreducible(Polynomial([1, 1]), 2) == false

	p = Polynomial([1, 1, 2, 0, 1, 2, 1])
	q = Polynomial([2, 1, 1])
	a = Polynomial([0, 1, 0, 1, 0, 0, 1])
	b = Polynomial([1, 0, 1])

	@test mod(rem(p, q), 3) == Polynomial([1, 1])
	@test mod(rem(a, b), 2) == Polynomial([1])

    @test isperfectpower(36) == true
	@test isprimepower(36) == false
	@test isperfectpower(9) == true
	@test isprimepower(9) == true
	@test isperfectpower(5) == false
	@test isprimepower(2) == true
	@test isprimepower(18) == false

	@test rref([1 0 1 1 1; 1 1 1 0 1; 0 1 1 1 1], 2) == [1 0 0 1 0; 0 1 0 1 0; 0 0 1 0 1]
	@test rref([1 0 1 0 1 0; 0 1 0 0 1 0; 1 1 1 1 1 1], 2) == [1 0 1 0 1 0; 0 1 0 0 1 0; 0 0 0 1 1 1]
	@test rref([1 1 0 2 3 1; 2 0 1 3 4 1; 1 2 2 1 4 3], 5, colswap=false) == [1 0 3 0 2 2; 0 1 2 0 1 1; 0 0 0 1 0 4]
	@test rref([1 1 0 2 3 1; 2 0 1 3 4 1; 1 2 2 1 4 3], 5, colswap=true) == [1 0 0 3 2 2; 0 1 0 2 1 1; 0 0 1 0 0 4]
	@test rref([1 2 0 1 2 1 2;  2 2 2 0 1 1 1; 1 0 1 1 2 1 2; 0 1 0 1 1 2 2], 3) == [1 0 0 0 2 2 2; 0 1 0 0 2 0 1; 0 0 1 0 1 0 2; 0 0 0 1 2 2 1]
	@test rref([0 0 0 0 0; 1 0 1 0 1; 0 1 0 1 1; 1 1 1 1 0], 2) == [1 0 1 0 1; 0 1 0 1 1; 0 0 0 0 0; 0 0 0 0 0]
	@test rref([1 1 1 0; 1 1 0 1; 0 0 1 1], 2) == [1 1 0 1; 0 0 1 1; 0 0 0 0]

	@test multiplication_table(2, 3) == Polynomial[Polynomial([0]) Polynomial([0]) Polynomial([0]) Polynomial([0]) Polynomial([0]) Polynomial([0]) Polynomial([0]) Polynomial([0]) Polynomial([0]); Polynomial([0]) Polynomial([1]) Polynomial([2]) Polynomial([0, 1]) Polynomial([1, 1]) Polynomial([2, 1]) Polynomial([0, 2]) Polynomial([1, 2]) Polynomial([2, 2]); Polynomial([0]) Polynomial([2]) Polynomial([1]) Polynomial([0, 2]) Polynomial([2, 2]) Polynomial([1, 2]) Polynomial([0, 1]) Polynomial([2, 1]) Polynomial([1, 1]); Polynomial([0]) Polynomial([0, 1]) Polynomial([0, 2]) Polynomial([0, 0, 1]) Polynomial([0, 1, 1]) Polynomial([0, 2, 1]) Polynomial([0, 0, 2]) Polynomial([0, 1, 2]) Polynomial([0, 2, 2]); Polynomial([0]) Polynomial([1, 1]) Polynomial([2, 2]) Polynomial([0, 1, 1]) Polynomial([1, 2, 1]) Polynomial([2, 0, 1]) Polynomial([0, 2, 2]) Polynomial([1, 0, 2]) Polynomial([2, 1, 2]); Polynomial([0]) Polynomial([2, 1]) Polynomial([1, 2]) Polynomial([0, 2, 1]) Polynomial([2, 0, 1]) Polynomial([1, 1, 1]) Polynomial([0, 1, 2]) Polynomial([2, 2, 2]) Polynomial([1, 0, 2]); Polynomial([0]) Polynomial([0, 2]) Polynomial([0, 1]) Polynomial([0, 0, 2]) Polynomial([0, 2, 2]) Polynomial([0, 1, 2]) Polynomial([0, 0, 1]) Polynomial([0, 2, 1]) Polynomial([0, 1, 1]); Polynomial([0]) Polynomial([1, 2]) Polynomial([2, 1]) Polynomial([0, 1, 2]) Polynomial([1, 0, 2]) Polynomial([2, 2, 2]) Polynomial([0, 2, 1]) Polynomial([1, 1, 1]) Polynomial([2, 0, 1]); Polynomial([0]) Polynomial([2, 2]) Polynomial([1, 1]) Polynomial([0, 2, 2]) Polynomial([2, 1, 2]) Polynomial([1, 0, 2]) Polynomial([0, 1, 1]) Polynomial([2, 0, 1]) Polynomial([1, 2, 1])]

	@test aresamearrays(list_span([2, 1, 1], [1, 1, 1], 3), [[0, 0, 0], [1, 1, 1], [2, 2, 2], [2, 1, 1], [0, 2, 2], [1, 0, 0], [1, 2, 2], [2, 0, 0], [0, 1, 1]])

	@test islinear([[0,0,0],[1,1,1],[1,0,1],[1,1,0]], 2) == false
	@test islinear([[0,0,0],[1,1,1],[1,0,1],[0,1,0]], 2) == true

	@test code_distance([[0,0,0,0,0],[1,0,1,0,1],[0,1,0,1,0],[1,1,1,1,1]]) == 2
	@test code_distance([[0,0,0,0,0],[1,1,1,0,0],[0,0,0,1,1],[1,1,1,1,1],[1,0,0,1,1],[0,1,1,0,0]]) == 1

	w1 = Word((1, 5, 2, 6)); w2 = Word([4, 2, 6, 7, 3, 10, 10]); w3 = Word(1, 3, 6, 2, 4)
	@test all(isword.([w1, w2, w3])) == true
	@test all(isabstractword.([w1, w2, w3])) == true
	@test isabstractword(1, 2, 3) == true
	@test all(isabstractword.([(1, 2, 3), [1, 2, 3]])) == true
	@test any(isword.([(1, 2, 3), [4, 5, 6]])) == false
	@test isword(1, 2, 3) == false
	@test aresamearrays(Alphabet("123"), deepsym([1, 2, 3]))
	@test aresamearrays(Alphabet([1, 2, 3]), deepsym([1, 2, 3]))
	@test aresamearrays(Alphabet(["1", "2", "3"]), deepsym([1, 2, 3]))
	# TODO: write test for CodeUniverse struct
	@test aresamearrays([i for i in CodeUniverseIterator(["a", "b", "c"], 3)], get_all_words(["a", "b", "c"], 3)) # implicitly tests CodeUniverseIterator
	@test aresamearrays(collect(CodeUniverseIterator(["a", "b", "c"], 4)), Tuple[(:a, :a, :a, :a), (:b, :a, :a, :a), (:c, :a, :a, :a), (:a, :b, :a, :a), (:b, :b, :a, :a), (:c, :b, :a, :a), (:a, :c, :a, :a), (:b, :c, :a, :a), (:c, :c, :a, :a), (:a, :a, :b, :a), (:b, :a, :b, :a), (:c, :a, :b, :a), (:a, :b, :b, :a), (:b, :b, :b, :a), (:c, :b, :b, :a), (:a, :c, :b, :a), (:b, :c, :b, :a), (:c, :c, :b, :a), (:a, :a, :c, :a), (:b, :a, :c, :a), (:c, :a, :c, :a), (:a, :b, :c, :a), (:b, :b, :c, :a), (:c, :b, :c, :a), (:a, :c, :c, :a), (:b, :c, :c, :a), (:c, :c, :c, :a), (:a, :a, :a, :b), (:b, :a, :a, :b), (:c, :a, :a, :b), (:a, :b, :a, :b), (:b, :b, :a, :b), (:c, :b, :a, :b), (:a, :c, :a, :b), (:b, :c, :a, :b), (:c, :c, :a, :b), (:a, :a, :b, :b), (:b, :a, :b, :b), (:c, :a, :b, :b), (:a, :b, :b, :b), (:b, :b, :b, :b), (:c, :b, :b, :b), (:a, :c, :b, :b), (:b, :c, :b, :b), (:c, :c, :b, :b), (:a, :a, :c, :b), (:b, :a, :c, :b), (:c, :a, :c, :b), (:a, :b, :c, :b), (:b, :b, :c, :b), (:c, :b, :c, :b), (:a, :c, :c, :b), (:b, :c, :c, :b), (:c, :c, :c, :b), (:a, :a, :a, :c), (:b, :a, :a, :c), (:c, :a, :a, :c), (:a, :b, :a, :c), (:b, :b, :a, :c), (:c, :b, :a, :c), (:a, :c, :a, :c), (:b, :c, :a, :c), (:c, :c, :a, :c), (:a, :a, :b, :c), (:b, :a, :b, :c), (:c, :a, :b, :c), (:a, :b, :b, :c), (:b, :b, :b, :c), (:c, :b, :b, :c), (:a, :c, :b, :c), (:b, :c, :b, :c), (:c, :c, :b, :c), (:a, :a, :c, :c), (:b, :a, :c, :c), (:c, :a, :c, :c), (:a, :b, :c, :c), (:b, :b, :c, :c), (:c, :b, :c, :c), (:a, :c, :c, :c), (:b, :c, :c, :c), (:c, :c, :c, :c)])

	@test sphere_covering_bound(5,7,3) == 215
	@test sphere_packing_bound(5,7,3) == 2693

	@test aresamearrays(construct_ham_matrix(3,2), [0 0 0 1 1 1 1; 0 1 1 0 0 1 1; 1 0 1 0 1 0 1])
	@test aresamearrays(construct_ham_matrix(3,3), [0 0 0 0 0 0 0 0 1 1 1 1 1; 0 0 1 1 1 2 2 2 0 0 0 1 1; 1 2 0 1 2 0 1 2 0 1 2 0 1])

	@test isperfect(11, 6, 5, 3) == true
	@test isperfect(23, 12, 7, 2) == true
	@test isperfect(23, 12, 7, 3) == false
	@test isperfect(11, 6, 5, 4) == false
	@test isgolayperfect(11, 6, 5, 3) == true
	@test isgolayperfect(23, 12, 7, 2) == true
	@test isgolayperfect(23, 12, 7, 3) == false
	@test isgolayperfect(11, 6, 5, 4) == false

	@test length(get_codewords(5, 5, 3)) ∈ [74:74...] # implicitly tests Word struct, and `get_codewords_random`
	@test length(get_codewords(4, 7, 3; m = 1)) ∈ [256:308...]
	@test length(get_codewords_greedy(5, 5, 3)) == 74
	randq, randn = rand(1:8, 2)
	@test length(get_all_words(randq, randn)) == big(randq)^randn
	@test aresamearrays(get_codewords([1 0 1 0; 0 1 1 1], 2), [[0, 0, 0, 0], [1, 0, 1, 0], [0, 1, 1, 1], [1, 1, 0, 1]])
	@test aresamearrays(get_codewords([1 0 0 1 1 0; 0 1 0 1 0 1; 0 0 1 0 1 1], 2), [[0, 0, 0, 0, 0, 0], [1, 0, 0, 1, 1, 0], [0, 1, 0, 1, 0, 1], [1, 1, 0, 0, 1, 1], [0, 0, 1, 0, 1, 1], [1, 0, 1, 1, 0, 1], [0, 1, 1, 1, 1, 0], [1, 1, 1, 0, 0, 0]])

	# note that the following tests are not to be checked using `aresamearrays` because their order matters
	@test syndrome([0, 2, 1, 2, 0, 1, 0], transpose(parity_check([1 0 0 0 2 2 2; 0 1 0 0 2 0 1; 0 0 1 0 1 0 2; 0 0 0 1 2 2 1], 3)), 3) == [0 0 0]
	@test parity_check([1 0 0 0 2 2 2; 0 1 0 0 2 0 1; 0 0 1 0 1 0 2; 0 0 0 1 2 2 1], 3) == [1 1 2 1 1 0 0; 1 0 0 1 0 1 0; 1 2 1 2 0 0 1]
	@test normal_form([1 2 0 1 2 1 2; 2 2 2 0 1 1 1; 1 0 1 1 2 1 2; 0 1 0 1 1 2 2], 3) == [1 0 0 0 2 2 2; 0 1 0 0 2 0 1; 0 0 1 0 1 0 2; 0 0 0 1 2 2 1]
	@test equivalent_code([1 2 0 1 2 1 2; 2 2 2 0 1 1 1; 1 0 1 1 2 1 2; 0 1 0 1 1 2 2], 3) == [1 0 0 0 2 2 2; 0 1 0 0 2 0 1; 0 0 1 0 1 0 2; 0 0 0 1 2 2 1]
	@test parity_check(normal_form([1 2 0 1 2 1 2; 2 2 2 0 1 1 1; 1 0 1 1 2 1 2; 0 1 0 1 1 2 2], 3), 3) == [1 1 2 1 1 0 0; 1 0 0 1 0 1 0; 1 2 1 2 0 0 1]
	@test isincode([0, 2, 1, 2, 0, 1, 0], transpose(parity_check([1 0 0 0 2 2 2; 0 1 0 0 2 0 1; 0 0 1 0 1 0 2; 0 0 0 1 2 2 1], 3)), 3) == true
	@test isincode([1, 0, 2, 2, 1, 2, 1], transpose(parity_check([1 0 0 0 2 2 2; 0 1 0 0 2 0 1; 0 0 1 0 1 0 2; 0 0 0 1 2 2 1], 3)), 3) == false
end # end runtests


text:
using Unitful, UnitfulAstro

import PhysicalConstants.CODATA2018: g_n

import ..SfUnits: Angle, to_angle, Acceleration

export kinetic_energy, distance, duration, acceleration, projectile_displacement, projectile_velocity,
	projectile_flight_time, projectile_apex, projectile_range, projectile_angle, projectile_angle_planetary

# TODO: rename these when the underlying compilation issues have been resolved
export projectile_range_planetary
# TODO: better names?
export projectile_peak_displacement_planetary

"""
	kinetic_energy(m::Unitful.Mass, v::Unitful.Velocity)

Compute the kinetic energy of a body with mass `m` travelling at velocity `v`.

No relativistic corrections are applied. Use `relativistic_kinetic_energy` if those are required.
"""
kinetic_energy(m::Unitful.Mass, v::Unitful.Velocity) = 0.5m*v^2 |> u"J"

"""
    distance(a::Unitful.Acceleration, t::Unitful.Time, initial_v::Unitful.Velocity = 0u"m/s")

Compute the distance travelled in time `t` by a body with initial velocity `initial_v` and uniform acceleration `a`.
"""
distance(a::Unitful.Acceleration, t::Unitful.Time, initial_v::Unitful.Velocity = 0u"m/s") = t*initial_v + a*t^2 / 2

"""
    duration(a::Unitful.Acceleration, d::Unitful.Length, initial_v::Unitful.Velocity = 0u"m/s")

Compute the time taken for a body with initial velocity `initial_v` and uniform acceleration `a` to travel a distance of `d`.
"""
function duration(a::Unitful.Acceleration, d::Unitful.Length, initial_v::Unitful.Velocity = 0u"m/s")
    if d < 0u"m"
		throw(DomainError(d, "Negative distances not allowed."))
	end

	if initial_v == 0u"m/s"
		return sqrt(2d / a)
	elseif a == 0u"m/s^2"
		if initial_v < 0u"m/s"
			throw(DomainError(d, "Negative initial velocities not allowed with zero acceleration."))
		end

		return d / initial_v
	end

    # at^2 + 2v_0t - 2d = 0
    # ax^2 + bx + c

    b = 2 * initial_v
    c = -2d

    discr = b^2 - 4*a*c

    t1 = (-b + sqrt(discr))/(2a)
    t2 = (-b - sqrt(discr))/(2a)

    # the negative results aren't wrong, per se, but as the parameters are
    # vector quantities the negative results aren't particularly interesting.

    if t1 < 0u"s"
        return t2
    elseif t2 < 0u"s"
        return t1
    else
        return min(t1, t2)
    end
end

"""
    acceleration(d::Unitful.Length, v_final::Unitful.Velocity, v_initial::Unitful.Velocity = 0u"m/s")

Given a distance `d`, initial velocity `v_initial and desired velocity `v_final` compute the required uniform acceleration
"""
acceleration(d::Unitful.Length, v_final::Unitful.Velocity, v_initial::Unitful.Velocity = 0u"m/s") = (v_final^2 - v_initial^2) / 2d

"""
    projectile_displacement(v_0::Unitful.Velocity, θ::Angle, t::Unitful.Time; g::Acceleration = g_n)

Horizontal and vertical displacement at time `t` of a projectile with initial velocity `v_0` and launch angle `θ` in a unitform gravitational acceleration `g`.
"""
projectile_displacement(v_0::Unitful.Velocity, θ::Angle, t::Unitful.Time; g::Acceleration = g_n) =
    [ v_0 * t * cos(θ), v_0 * t * sin(θ) - 0.5g * t^2 ]

"""
    projectile_displacement(v_0::Unitful.Velocity, t::Unitful.Time; g::Acceleration = g_n)

Vertical displacement at time `t` of a projectile launched vertically with initial velocity `v_0` in a unitform gravitational acceleration `g`.
"""
projectile_displacement(v_0::Unitful.Velocity, t::Unitful.Time; g::Acceleration = g_n) = projectile_displacement(v_0, 90u"°", t, g = g)[2]

"""
    projectile_velocity(v_0::Unitful.Velocity, θ::Angle, t::Unitful.Time; g::Acceleration = g_n)

Horizontal and vertical velocity at time `t` of a projectile with initial velocity `v_0` and launch angle `θ` in a unitform gravitational acceleration `g`.
"""
projectile_velocity(v_0::Unitful.Velocity, θ::Angle, t::Unitful.Time; g::Acceleration = g_n) =
    [ v_0 * cos(θ), v_0 * sin(θ) - g * t ]

"""
    projectile_velocity(v_0::Unitful.Velocity, t::Unitful.Time; g::Acceleration = g_n)

Velocity at time `t` of a projectile launched vertically with initial velocity `v_0` in a unitform gravitational acceleration `g`.
"""
projectile_velocity(v_0::Unitful.Velocity, t::Unitful.Time; g::Acceleration = g_n) = projectile_velocity(v_0, 90u"°", g, t)[2]

"""
    projectile_flight_time(v_0::Unitful.Velocity, θ::Angle = 90u"°"; g::Acceleration = g_n)

Compute total flight time for a projectile launched at angle `θ` with initial velocity `v_0` in a uniform gravitation acceleration `g`.

Flight is considered to be complete when vertical displacement returns to zero.
"""
projectile_flight_time(v_0::Unitful.Velocity, θ::Angle = 90u"°"; g::Acceleration = g_n) = 2v_0 * sin(θ) / g |> u"s"

"""
    projectile_apex(v_0::Unitful.Velocity, θ::Angle = 90u"°"; g::Acceleration = g_n)

Horizontal and vertical displacement of a projectile launched at angle `θ` with initial velocity `v_0` in a uniform gravitation acceleration `g` at the apex of its trajectory.

Default launch angle is 90° for maximum height.
"""
projectile_apex(v_0::Unitful.Velocity, θ::Angle = 90u"°"; g::Acceleration = g_n) = [ v_0^2 * sin(2θ) / 2g,  v_0^2 * sin(θ)^2 / 2g ]

"""
    projectile_range(v_0::Unitful.Velocity, θ::Angle = 45u"°"; g::Acceleration = g_n)

Compute range of a projectile launched at angle `θ` with initial velocity `v_0` in a uniform gravitation acceleration `g`.

Default launch angle is 45° for maximum range.
"""
projectile_range(v_0::Unitful.Velocity, θ::Angle = 45u"°"; g::Acceleration = g_n) = (v_0^2 * sin(2θ)) / g

"""
    projectile_range(v_0::Unitful.Velocity, y_0::Unitful.Length, θ::Angle = 45u"°"; g::Acceleration = g_n)

Compute range of a projectile launched from altitude `y_0` at angle `θ` with initial velocity `v_0` in a uniform gravitation acceleration `g`.

Range is reached when vertical displacement is zero. Default launch angle is 45° for maximum range.
"""
projectile_range(v_0::Unitful.Velocity, y_0::Unitful.Length, θ::Angle = 45u"°"; g::Acceleration = g_n) =
	(v_0 * cos(θ) / g) * (v_0 * sin(θ) + sqrt((v_0 * sin(θ))^2 + 2g * y_0))

"""
    projectile_range_planetary(v_0::Unitful.Velocity, r_planet::Unitful.Length = 1u"Rearth", θ::Angle = 45u"°"; g::Acceleration = g_n)

Compute range of a projectile launched from altitude `y_0` at angle `θ` with initial velocity `v_0` on a spherical planet with radius `r_planet` and surface gravity `g`.

NOTE: this should be called `projectile_range`, but weird internal compiler errors prevent that with this version of Julia (1.5.3).
"""
function projectile_range_planetary(v_0::Unitful.Velocity, θ::Angle = 45u"°", r_planet::Unitful.Length = 1u"Rearth"; g::Acceleration = g_n)
	v_rat2 = (v_0 / sqrt(r_planet * g))^2

	if v_rat2 > 1
		throw(DomainError("Initial velocity exceeds orbital velocity; Range undefined"))
	end

	a = v_0^2 * sin(2θ) / g
	b = sqrt(1 - (2 - v_rat2) * v_rat2 * cos(θ)^2)

	return a / b |> u"m"
end

function projectile_peak_displacement_planetary(v_0::Unitful.Velocity, θ::Angle = 45u"°", r_planet::Unitful.Length = 1u"Rearth"; g::Acceleration = g_n)
	v_rat2 = (v_0 / sqrt(r_planet * g))^2

	if v_rat2 > 1
		throw(DomainError("Initial velocity exceeds orbital velocity; Range undefined"))
	end

	a = v_0^2 * sin(θ) / g
	b = 1 - v_rat2 + sqrt(1 - (2 - v_rat2) * v_rat2 * cos(θ)^2)

	return a / b |> u"m"
end

"""
    projectile_angle(v_0::Unitful.Velocity, d::Unitful.Length; g::Acceleration = g_n)

Compute the possible launch angles for a projectile to reach a horizontal range of `d` given initial velocity `v_0` in a uniform gravitation acceleration `g`.

The shallow angle is the first result, the steep angle the second. A domain error is raised if the target is out of range.
"""
function projectile_angle(v_0::Unitful.Velocity, d::Unitful.Length; g::Acceleration = g_n)
	k = (g * d) / v_0^2

	if (k > 1)
		throw(DomainError("Target out of range; no solutions"))
	end

	return ( 0.5asin(k) * 1u"rad" |> u"°", 45u"°" + (0.5acos(k) * 1u"rad" |> u"°") )
end

"""
    projectile_angle(v_0::Unitful.Velocity, d::Unitful.Length, y::Unitful.Length; g::Acceleration = g_n)

Compute the possible launch angles for a projectile to reach a horizontal range of `d` and altitude of `y` relative to the starting point given initial velocity `v_0` in a uniform gravitation acceleration `g`.

The shallow angle is the first result, the steep angle the second. A domain error is raised if the target is out of range.
"""
function projectile_angle(v_0::Unitful.Velocity, d::Unitful.Length, y::Unitful.Length; g::Acceleration = g_n)
	det = v_0^4 - g * (g * d^2 + 2y * v_0^2)

	if (ustrip(det) < 0)
		throw(DomainError("Target out of range; no solutions"))
	end

	a = atan((v_0^2 + sqrt(det)) / (g * d)) * 1u"rad" |> u"°"
	b = atan((v_0^2 - sqrt(det)) / (g * d)) * 1u"rad" |> u"°"

	if a < b
		return (a, b)
	else
		return (b, a)
	end
end

"""
    projectile_optimum_angle(v_0::Unitful.Velocity, r_planet::Unitful.Length = 1u"Rearth"; g::Acceleration = g_n)
"""
function projectile_optimum_angle(v_0::Unitful.Velocity, r_planet::Unitful.Length = 1u"Rearth"; g::Acceleration = g_n)
	v_rat2 = (v_0 / sqrt(r_planet * g))^2

	if v_rat2 > 0
		throw(DomainError("Initial velocity exceeds orbital velocity; optimum angle undefined"))
	end

	return 0.5acos(v_rat2 / (2 - v_rat2)) * u"rad" |> u"°"
end

ext: jl

max_stars_repo_name: icantbelieveitsnotphysics/SfPhysics

max_issues_repo_name: icantbelieveitsnotphysics/SfPhysics

max_forks_repo_name: icantbelieveitsnotphysics/SfPhysics


text:
#=

  Euler #49 in Julia.

  Problem 49:
  """
  The arithmetic sequence, 1487, 4817, 8147, in which each of the terms
  increases by 3330, is unusual in two ways: (i) each of the three terms are
  prime, and, (ii) each of the 4-digit numbers are permutations of one another.

  There are no arithmetic sequences made up of three 1-, 2-, or 3-digit primes,
  exhibiting this property, but there is one other 4-digit increasing sequence.

  What 12-digit number do you form by concatenating the three terms
  in this sequence?
  """

  This Julia program was created by Hakan Kjellerstrand, hakank@gmail.com
  See also my Julia page: http://www.hakank.org/julia/

=#

include("Euler.jl")

function check_perms(n, diff)
    allperms = all_permutations2(digits(n))
    if length(allperms) > 0
        p1 = get_element(n, allperms, diff)
        if p1 !== nothing
            p2 = get_element(p1, allperms, diff)
            if p2 !== nothing
                return [n, p1, p2]
            end
        end
    end

    return nothing
end

function get_element(n, ll, diff)
    for p in ll
        pp = parse(Int,join(p,""))
        # if isPrime(pp) && pp > n && pp-n == diff
        if pp > n && pp-n == diff && isPrime(pp)
            return pp
        end
    end
    return nothing
end

# 0.10001170s
function euler49a()
    diff = 3330
    res = 0
    for n in 1001:2:9999
        if n != 1487 && isPrime(n)
            c = check_perms(n, diff)
            if c !== nothing
                res = c
                break
            end
        end
    end

    return parse(Int,join(res,""))

end

run_euler(euler49a)

max_stars_repo_name: Wikunia/hakank

max_issues_repo_name: Wikunia/hakank

max_forks_repo_name: Wikunia/hakank


text:
function _precompile_()
    ccall(:jl_generating_output, Cint, ()) == 1 || return nothing
    isdefined(Core, Symbol("#Type##kw")) && precompile(Tuple{getfield(Core, Symbol("#Type##kw")), NamedTuple{(:libc, :compiler_abi), Tuple{Nothing, Pkg.BinaryPlatforms.CompilerABI}}, Type{Pkg.BinaryPlatforms.FreeBSD}, Symbol})
    isdefined(Core, Symbol("#Type##kw")) && precompile(Tuple{getfield(Core, Symbol("#Type##kw")), NamedTuple{(:libgfortran_version, :libstdcxx_version, :cxxstring_abi), Tuple{Nothing, Nothing, Symbol}}, Type{Pkg.BinaryPlatforms.CompilerABI}})
end


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

# Deterministic equivalent optimization #
# ===================================== #
function optimize!(structure::DeterministicEquivalent, optimizer::MOI.AbstractOptimizer, x₀::AbstractVector)
    # Sanity check
    backend(structure.model) === optimizer || error("Stochastic program optimizer has not been connected to the deterministically equivalent problem.")
    # Crash if supported
    if MOI.supports(optimizer, MOI.VariablePrimalStart(), MOI.VariableIndex)
        for (i, dvar) in enumerate(all_decision_variables(structure.model, 1))
            MOI.set(optimizer, MOI.VariablePrimalStart(), index(dvar), x₀[i])
        end
    end
    # Run standard MOI optimization procedure
    MOI.optimize!(optimizer)
    return nothing
end

function supports_structure(::MOI.AbstractOptimizer, ::DeterministicEquivalent)
    return true
end

function check_loadable(::MOI.AbstractOptimizer, ::DeterministicEquivalent)
    return nothing
end

function load_structure!(::MOI.AbstractOptimizer, ::DeterministicEquivalent, ::AbstractVector)
    return nothing
end

function set_master_optimizer!(structure::DeterministicEquivalent, optimizer)
    # Ensure decision bridges are added
    for bridge_type in structure.model.bridge_types
        JuMP._moi_add_bridge(structure.model.moi_backend, bridge_type)
    end
    return nothing
end

function set_master_optimizer_attribute!(structure::DeterministicEquivalent, attr::MOI.AbstractOptimizerAttribute, value)
    MOI.set(backend(structure.model), attr, value)
    return nothing
end

function set_subproblem_optimizer!(structure::DeterministicEquivalent, optimizer)
    return nothing
end

function set_subproblem_optimizer_attribute!(::DeterministicEquivalent, ::MOI.AbstractOptimizerAttribute, value)
    return nothing
end

function optimize!(structure::DeterministicEquivalent, optimizer::AbstractStructuredOptimizer)
    # Sanity check
    supports_structure(optimizer, typeof(structure)) || throw(UnsupportedStructure{typeof(optimizer), typeof(structure)}())
    # Ensure that no decisions are fixed in the first stage
    untake_decisions!(structure.model, structure.decision_variables[1])
    # Load structure
    load_structure!(optimizer, structure)
    # Run structure-exploiting optimization procedure
    MOI.optimize!(optimizer)
    return nothing
end

function cache_solution!(stochasticprogram::StochasticProgram{2}, structure::DeterministicEquivalent{2}, optimizer::MOI.AbstractOptimizer)
    cache = solutioncache(stochasticprogram)
    # Cache main solution
    variables = decision_variables_at_stage(stochasticprogram, 1)
    constraints = decision_constraints_at_stage(stochasticprogram, 1)
    cache[:solution] = SolutionCache(backend(structure.model), variables, constraints)
    # Cache first-stage solution
    cache[:node_solution_1] = SolutionCache(backend(structure.model), variables, constraints)
    try
        Q = MOIU.eval_variables(structure.decisions.stage_objectives[1][1][2]) do idx
            return MOI.get(backend(structure.model), MOI.VariablePrimal(), idx)
        end
        cache[:node_solution_1].modattr[MOI.ObjectiveValue()] = Q
        cache[:node_solution_1].modattr[MOI.DualObjectiveValue()] = Q
    catch
    end
    # Cache scenario-dependent solutions
    variables = decision_variables_at_stage(stochasticprogram, 2)
    constraints = decision_constraints_at_stage(stochasticprogram, 2)
    for scenario_index in 1:num_scenarios(stochasticprogram, 2)
        key = Symbol(:node_solution_2_, scenario_index)
        cache[key] = SolutionCache(backend(structure.model))
        # Model attributes are shared
        cache_model_attributes!(cache[key], backend(structure.model))
        # Variables/constraints are scenario-dependent
        cache_variable_attributes!(cache[key], structure, variables, 2, scenario_index)
        cache_constraint_attributes!(cache[key], structure, constraints, 2, scenario_index)
        # Cache subobjective
        if cache[key].modattr[MOI.TerminationStatus()] == MOI.OPTIMAL
            Q = 0.0
            try
                Qᵢ = MOIU.eval_variables(structure.decisions.stage_objectives[2][i][2]) do idx
                    return MOI.get(backend(structure.model), MOI.VariablePrimal(), idx)
                end
                cache[key].modattr[MOI.ObjectiveValue()] = Qᵢ
                cache[key].modattr[MOI.DualObjectiveValue()] = Qᵢ
                Q += Qᵢ
            catch
            end
            # Correct first-stage objective
            cache[:node_solution_1].modattr[MOI.ObjectiveValue()] -= Q
        end
    end
    return nothing
end


text:
# MIT License
# Copyright (c) 2017: Xavier Gandibleux, Anthony Przybylski, Gauthier Soleilhac, and contributors.
struct solution
    pb::problem #The problem for which it's a solution
    x::BitArray #Variable states
    obj_1::Int #Value on first objective
    obj_2::Int #Value on second objective
    weight::Int #Weight of the solution
end

#Constructor : empty solution
solution(p::problem) = solution(p, falses(size(p)), 0, 0, 0)

#Getters
obj_1(s::solution) = s.obj_1 + s.pb.min_profit_1
obj_2(s::solution) = s.obj_2 + s.pb.min_profit_2
obj(s::solution) = obj_1(s), obj_2(s)
weight(s::solution) = s.weight + s.pb.ω

full_variable_vector(s::solution) = begin
    res = zeros(Int,s.pb.psize)
    for i in s.pb.C1
        res[i] = 1
    end
    for i = 1:size(s.pb)
        if s.x[i] == true
            res[s.pb.variables[i]] = 1
        end
    end
    return res
end

print_vars(s::solution) = String([i==1 ? '1' : '0' for i in full_variable_vector(s)])

#Print function
Base.show(io::IO, s::solution) =
    # print(io, "sol(p1=$(obj_1(s)), p2=$(obj_2(s)), w=$(weight(s)))")
    print(io, "sol(p1=",obj_1(s),", p2=",obj_2(s),") w=",weight(s)," : " ,print_vars(s))

#Helper function to calculate an objective value or an accumulated weight
function dot(v::Vector{Int}, status::BitArray, indices::AbstractArray{Int,1})
    #@assert length(v) >= length(status) == length(indices)
    res = 0
    for i = 1:length(status)
        @inbounds res += status[i] * v[indices[i]]
    end
    res
end

#Relation operators
dominates(a1::Int,a2::Int,b1::Int,b2::Int) = (a1 > b1 && a2 >= b2) || (a1 >= b1 && a2 > b2)
dominates(yr::Tuple{Int,Int}, ys::Tuple{Int,Int}) = dominates(yr[1], yr[2], ys[1], ys[2])

ideal(yr::Tuple{Int,Int}, ys::Tuple{Int,Int}) = max(yr[1],ys[1]), max(yr[2],ys[2])
nadir(yr::Tuple{Int,Int}, ys::Tuple{Int,Int}) = min(yr[1],ys[1]), min(yr[2],ys[2])

import Base.==
==(a::solution, b::solution) = a.obj_1==b.obj_1 && a.obj_2==b.obj_2 && a.weight==b.weight && a.x==b.x
import Base.<
import Base.<=
<=(a::solution,b::solution) = dominates(obj(b),obj(a))
<=(a, b::solution) = dominates(obj(b), a)
<=(a::solution, b) = dominates(b, obj(a))

ideal(a::solution,b::solution) = max(obj_1(a), obj_1(b)) , max(obj_2(a), obj_2(b))
nadir(a::solution,b::solution) = nadir(obj(a), obj(b))

max_stars_repo_name: vOptSolver/vOptSpecific

max_issues_repo_name: vOptSolver/vOptSpecific

max_forks_repo_name: vOptSolver/vOptSpecific


text:
@testset "Paired read datastores" begin
    function get_fastq_seqs(::Type{A}, file) where {A<:DNAAlphabet}
        seqs = map(open(FASTQ.Reader, file) do rdr
            collect(rdr)
        end) do rec
            FASTQ.sequence(LongSequence{A}, rec)
        end
        for seq in seqs
            if length(seq) > 300
                resize!(seq, 300)
            end
        end
        return seqs
    end

    function check_round_trip(::Type{A}, R1, R2) where {A<:DNAAlphabet}
        r1_seqs = get_fastq_seqs(A, R1)
        r2_seqs = get_fastq_seqs(A, R2)

        fqa = open(FASTQ.Reader, R1)
        fqb = open(FASTQ.Reader, R2)

        ds = PairedReads{A}(fqa, fqb, "ecoli-pe", "ecoli-pe", 250, 300, 0, FwRv)

        ds_seqs = collect(ds)

        return ds_seqs[1:2:end] == r1_seqs && ds_seqs[2:2:end] == r2_seqs
    end

    function check_show(ds, msg)
        buf = IOBuffer()
        show(buf, ds)
        return String(take!(buf)) == msg
    end

    function run_checks(::Type{A}) where {A<:DNAAlphabet}
        @test check_round_trip(A, "ecoli_tester_R1.fastq", "ecoli_tester_R2.fastq")
        ds = open(PairedReads{A}, "ecoli-pe.prseq")
        @test ReadDatastores.deduce_datastore_type("ecoli-pe.prseq") == PairedReads{A}
        @test ReadDatastores.name(ds) == Symbol("ecoli-pe")
        @test ReadDatastores.max_read_length(ds) == 300
        @test ReadDatastores.orientation(ds) == FwRv
        @test check_show(ds, "Paired Read Datastore 'ecoli-pe': 20 reads (10 pairs)")
        @test firstindex(ds) == 1
        @test lastindex(ds) == 20
        @test Base.IteratorSize(ds) == Base.HasLength()
        @test Base.IteratorEltype(ds) == Base.HasEltype()
        @test Base.eltype(ds) == LongSequence{A}
        @test_throws BoundsError ds[100]
        @test_throws BoundsError buffer(ds)[100]
        @test_throws ErrorException buffer(ds, 1)
        @test_throws ErrorException buffer(ds, ReadDatastores._bytes_per_read(ds) - 1)
        @test_throws BoundsError load_sequence!(ds, 100, LongSequence{A}())

        @test collect(ds) == collect(buffer(ds)) == open(PairedReads{A}, "ecoli-pe.prseq") do ds
            collect(ds)
        end
        @test ds[5] == buffer(ds)[5] == load_sequence!(ds, 5, LongSequence{A}()) == load_sequence!(buffer(ds), 5, LongSequence{A}())
    end

    run_checks(DNAAlphabet{4})
    run_checks(DNAAlphabet{2})

    @test_throws ReadDatastores.DatastoreVersionError{PairedReads{DNAAlphabet{2}}} open(PairedReads{DNAAlphabet{2}}, "ecoli-paired-old.prseq")
    @test_throws ReadDatastores.DatastoreEncodingError{PairedReads{DNAAlphabet{4}}} open(PairedReads{DNAAlphabet{4}}, "ecoli-pe.prseq")
end


text:
function CursorSchema(stmt::QueryStmt)

    columns_count = ncol(stmt)
    column_query_info = columns_info(stmt)

    column_names_index = Dict{String, Int}()

    for (column_index, query_info) in enumerate(column_query_info)
        column_names_index[column_name(query_info)] = column_index
    end

    return CursorSchema(column_query_info, column_names_index)
end

@inline ncol(schema::CursorSchema) = length(schema.column_names_index)
@inline ncol(row::ResultSetRow) = ncol(row.cursor)
@inline ncol(cursor::Cursor) = ncol(cursor.stmt)
@inline ncol(rs::ResultSet) = ncol(rs.schema)

@inline nrow(rs::ResultSet) = length(rs.rows)

@inline Base.isempty(rs::ResultSet) = isempty(rs.rows)
@inline Base.size(rs::ResultSet) = ( nrow(rs), ncol(rs) )

@inline stmt(cursor::Cursor) = cursor.stmt

function Cursor(stmt::QueryStmt)
    schema = CursorSchema(stmt)
    return Cursor(stmt, schema)
end

@static if VERSION < v"0.7-"
    # Iteration protocol for Julia v0.6

    #=
    for i = I   # or  "for i in I"
        # body
    end
    is translated into:

    state = start(I)
    while !done(I, state)
        (i, state) = next(I, state)
        # body
    end
    =#

    function Base.start(cursor::Cursor) :: FetchResult
        return fetch(cursor.stmt)
    end

    function Base.done(cursor::Cursor, state::FetchResult)
        return !state.found
    end

    function Base.next(cursor::Cursor, state::FetchResult)
        current_row = ResultSetRow(cursor)
        next_state = fetch(cursor.stmt)
        return (current_row, next_state)
    end

else
    # Iteration protocol for Julia v0.7 and v1.0

    #=
    for i in iter   # or  "for i = iter"
        # body
    end

    is translated into:

    next = iterate(iter)
    while next !== nothing
        (i, state) = next
        # body
        next = iterate(iter, state)
    end
    =#

    function Base.iterate(cursor::Cursor, nil::Nothing=nothing)
        fetch_result = fetch(cursor.stmt)
        if fetch_result.found
            current_row = ResultSetRow(cursor)
            return (current_row, nothing)
        else
            return nothing
        end
    end
end

function query(f::Function, conn::Connection, sql::AbstractString;
               scrollable::Bool=false,
               tag::AbstractString="",
               fetch_array_size::Integer=ORA_DEFAULT_FETCH_ARRAY_SIZE,
               exec_mode::OraExecMode=ORA_MODE_EXEC_DEFAULT)

    stmt(conn, sql, scrollable=scrollable, tag=tag, fetch_array_size=fetch_array_size) do stmt
        execute(stmt, exec_mode=exec_mode)
        cursor = Cursor(stmt)
        f(cursor)
    end
    nothing
end

function execute_and_fetch_all!(stmt::Stmt; exec_mode::OraExecMode=ORA_MODE_EXEC_DEFAULT) :: ResultSet
    local schema::CursorSchema
    rows = Vector{ResultSetRow}()

    query(stmt, exec_mode=exec_mode) do cursor
        schema = cursor.schema

        for row in cursor
            push!(rows, row)
        end
    end

    return ResultSet(schema, rows)
end

function query(f::Function, stmt::Stmt; exec_mode::OraExecMode=ORA_MODE_EXEC_DEFAULT)
    execute(stmt, exec_mode=exec_mode)
    f(Cursor(stmt))
    nothing
end

function query(stmt::Stmt; exec_mode::OraExecMode=ORA_MODE_EXEC_DEFAULT) :: ResultSet
    return execute_and_fetch_all!(stmt, exec_mode=exec_mode)
end

function query(conn::Connection, sql::AbstractString;
               scrollable::Bool=false,
               tag::AbstractString="",
               fetch_array_size::Integer=ORA_DEFAULT_FETCH_ARRAY_SIZE,
               exec_mode::OraExecMode=ORA_MODE_EXEC_DEFAULT) :: ResultSet

    local rs::ResultSet

    stmt(conn, sql, scrollable=scrollable, tag=tag, fetch_array_size=fetch_array_size) do stmt
        rs = execute_and_fetch_all!(stmt, exec_mode=exec_mode)
    end

    return rs
end

function ResultSetRow(cursor::Cursor)
    data = Vector{Any}()

    # parse data from Cursor
    for column_index in 1:ncol(cursor)
        oracle_value = query_oracle_value(stmt(cursor), column_index)
        push!(data, parse_oracle_value(oracle_value))
    end

    return ResultSetRow(cursor.schema, data)
end

Base.length(row::ResultSetRow) = length(row.data)
Base.isempty(row::ResultSetRow) = isempty(row.data)

@inline function Base.getindex(row::ResultSetRow, column_index::Integer)
    check_inbounds(row, column_index)
    @inbounds return row.data[column_index]
end

@inline function Base.getindex(row::ResultSetRow, column_name::AbstractString)
    check_inbounds(row, column_name)
    column_index = row.schema.column_names_index[column_name]
    @inbounds return row.data[column_index]
end

@inline function check_inbounds(schema::CursorSchema, column::Integer)
    @assert 0 < column <= ncol(schema) "Column $column not found."
end

@inline function check_inbounds(schema::CursorSchema, column::AbstractString)
    @assert haskey(schema.column_names_index, column) "Column $column not found."
end

@inline check_inbounds(row::ResultSetRow, column::Union{AbstractString,Integer}) = check_inbounds(row.schema, column)

@inline function check_inbounds(rs::ResultSet, row::Integer, column::Union{AbstractString, Integer})
    check_inbounds(rs.schema, column)
    @assert 0 < row <= nrow(rs) "Row $row not found."
end

has_possibly_more_rows(r::FetchRowsResult) = Bool(r.more_rows)

Base.IteratorSize(::Cursor) = Base.SizeUnknown()
Base.eltype(::Cursor) = ResultSetRow

@static if VERSION < v"0.7-"
    function Base.collect(cursor::Cursor)
        result = Vector{eltype(cursor)}()
        for row in cursor
            push!(result, row)
        end
        return result
    end
end

function fetchrow(stmt::QueryStmt) :: Union{Nothing, ResultSetRow}
    fetch_result = fetch(stmt)

    if fetch_result.found
        return ResultSetRow(Cursor(stmt))
    else
        return nothing
    end
end

function Base.getindex(rs::ResultSet, row::Integer, column::Union{AbstractString,Integer})
    check_inbounds(rs, row, column)
    return rs.rows[row][column]
end

function Base.getindex(rs::ResultSet, ::Colon, column::Union{AbstractString,Integer})
    column_data = Vector()

    !isempty(rs) && check_inbounds(rs.rows[1], column)

    for row in rs.rows
        push!(column_data, row[column])
    end

    return column_data
end

function _lastindex(rs::ResultSet, d::Integer)
    @assert d == 1 || d == 2 "Invalid dimension for ResultSet: $d."

    # rows
    if d == 1
        return nrow(rs)
    else
        # columns
        return ncol(rs)
    end
end

@static if VERSION < v"0.7-"
    Base.size(rs::ResultSet, d::Integer) = _lastindex(rs, d)
else
    Base.lastindex(rs::ResultSet, d::Integer) = _lastindex(rs, d)
end

function Base.getindex(rs::ResultSet, row_range::UnitRange{T}, column::Union{AbstractString,Integer}) where {T<:Integer}
    column_data = Vector()

    for r in row_range
        check_inbounds(rs, r, column)
        push!(column_data, rs[r, column])
    end

    return column_data
end


text:
# # Direct and inverse FFTW

# Here we test the methods, in FFTW, to be used in our pseudo-spectral code. It is mostly a sanity check.

# We start by loading the packages we need:

using FFTW
using WGLMakie
using Test

# ## The spatial domain and its discretization

# We consider a square domain of sides $L = 2\pi$, for which the smallest wavenumber is $\kappa_0 = 2\pi/L$.

L = 2π
κ₀ = 2π/L
nothing

# We set the number $N$ of points for the mesh in each direction, yielding a mesh $(x_i, y_j)_{i, j = 1,\ldots ,N}$, with $x_N = y_N = L$, and steps $x_{i+1} - x_i = L/N$, and $y_{j+1} - y_j = L/N$. Due to the periodicity, we don't need to store the values corresponding to $i = j = 0$.

N = 128
x = y = (L/N):(L/N):L
nothing

# We may visualize the grid with a scatter plot, although if the mesh is too thin, we won't quite see the details. If using GLMakie or WGLMakie, one can zoom in for a detailed view.

fig, ax, plt = scatter(vec(x .* one.(y)'), vec(one.(x) .* y'))

fig

# ## A vorticity function for testing

# In order to test the methods from [FFTW.jl](http://www.fftw.org), we define a certain vorticity function and its derivatives

vort = sin.(one.(y) * x') .* cos.(3y * one.(x)')
dx_vort = cos.(one.(y) * x') .* cos.(3y * one.(x)')
dy_vort = - 3 * sin.(one.(y) * x') .* sin.(3y * one.(x)')
dd_vort = - sin.(one.(y) * x') .* cos.(3y * one.(x)') - 9 * sin.(one.(y) * x') .* cos.(3y * one.(x)')

nothing

# We may visualize the vorticity as a surface:

fig = Figure(backgroundcolor = RGBf(0.98, 0.98, 0.98))

ax = Axis3(fig[1, 1], xlabel = "x", ylabel = "y", zlabel = "vorticity level", title = "Vorticity graph")

surface!(ax, x, y, vort, colormap = :berlin)

fig

# Or, better yet, as a heatmap:

fig, ax, plt = heatmap(x, y, vort, colormap = :berlin)

ax = Axis(fig[1, 1], xlabel = "x", ylabel = "y", title = "Vorticity heatmap")

fig

# ## Testing the direct and inverse discrete Fourier transform

# We check going back and forth with the real fft and inverse real fft:

@testset "check composition `irfft ∘ rfft`" begin
    vort_hat = rfft(vort)
    vortback = irfft(vort_hat, N)
    @test vort ≈ vortback
end
nothing

# We can also visualize the excited modes:

κ_x = κ₀ * (0:1/N:1/2)
κ_y = κ₀ * (0:1/N:1-1/N)

vort_hat = rfft(vort)

fig, ax, plt = heatmap(κ_x, κ_y, abs.(vort_hat).^2, colormap = :berlin)

ax = Axis(fig[1, 1], xlabel = "k_x", ylabel = "k_y", title = "Enstrophy spectrum heatmap")

fig

# We have excited modes $(\pm 1, \pm 3)$. Due to the reality condition of the vorticity, the Fourier spectrum has an Hermitian symmetry around the origin. Using the real Fourier transform, only half of the modes need to be stored. In this case, only modes $(1, \pm 3)$ are retained. Moreover, the negative modes are shifted above. Due to the one-indexing of Julia, this means that modes $(\pm a, \pm b)$ are represented by indices $[a + 1, b + 1]$ and $[a + 1, N + 1 - b]$. In our case, the excited modes are associated with

vort_hat[4, 2]

# and

vort_hat[4, 128]

# These are the only excited modes, as we can check:

for i in 1:div(N, 2) + 1, j in 1:N
    if abs(vort_hat[i, j])^2 > eps()
        println("vort_hat[$i, $j] = $(vort_hat[i, j])")
    end
end

# ## Differential operators

# In order to check the derivatives in spectral space, we define the following operators. Actually, they are just vectors, since derivatives in spectral space act as diagonal operators. Hence, a straighforward Hadamard product suffices.

Dx_hat = im * κ₀ * [ifelse(k1 ≤ N/2 + 1, k1 - 1, k1 - 1 - N) for k2 in 1:N/2+1, k1 in 1:N]
Dy_hat = im * κ₀ * [k2 - 1 for k2 in 1:N/2+1, k1 in 1:N]

Delta_hat = - κ₀^2 * [
    ifelse(k1 ≤ N/2 + 1, (k1 - 1)^2 + (k2 - 1)^2, (k1 - 1 - N)^2 + (k2 - 1)^2)
    for k2 in 1:N/2+1, k1 in 1:N
]

Hu_hat = - Dy_hat ./ Delta_hat
Hu_hat[1, 1] = 0.0
Hv_hat = Dx_hat ./ Delta_hat
Hv_hat[1, 1] = 0.0

nothing

# ## Testing the spectral operators

# Now we are ready to test the differentiation in spectral space, with the operators just defined.

@testset "Check operators" begin
    vort_hat = rfft(vort)
    u_hat = Hu_hat .* vort_hat
    v_hat = Hv_hat .* vort_hat
    ## check derivative dx in spectral space:
    @test irfft(Dx_hat .* vort_hat, N) ≈ dx_vort
    ## check derivative dy in spectral space:
    @test irfft(Dy_hat .* vort_hat, N) ≈ dy_vort
    ## check Laplacian in spectral space:
    @test Delta_hat ≈ Dx_hat.^2 .+ Dy_hat.^2
    @test irfft(Delta_hat .* vort_hat, N) ≈ dd_vort
    ## check recovering velocity field from vorticity
    @test Dx_hat .* v_hat - Dy_hat .* u_hat ≈ vort_hat
end
nothing

# ## One-mode steady state

# For forced-periodic flows, when the forcing function contains a single Fourier mode, there is a corresponding steady state vorticity function corresponding also to that single-mode in Fourier space, only with a different amplitude.

# In order to construct such pairs of forcing mode / steady state, we set the viscosity of the flow, choose the mode to be forced and define the strength of that force:

ν = 1.0e-0 # viscosity
κ = (x = 1, y = 2) # forced mode
α = (re = 0.1, im = 0.05) # strength

nothing

# With these parameters, we find the curl `g_steady` of the forcing term and the vorticity `vort_steady` of the corresponding steady-state:

g_steady = ν * (
    2α.re * κ₀^4 * sum(abs2, κ)^2 * cos.(κ₀ * (κ.x * one.(y) * x' + κ.y * y * one.(x)'))
    - 2α.im * κ₀^4 * sum(abs2, κ)^2 * sin.(κ₀ * (κ.x * one.(y) * x' + κ.y * y * one.(x)'))
)

vort_steady = 2κ₀^2 * sum(abs2, κ) * (
    α.re * cos.(κ₀ * (κ.x * one.(y) * x' + κ.y * y * one.(x)'))
    - α.im * sin.(κ₀ * (κ.x * one.(y) * x' + κ.y * y * one.(x)'))
)
nothing

# Here are some visualizations of the curl of the forcing term and of the vorticity

heatmap(x, y, g_steady, xlabel="x", ylabel="y", title="Curl of the forcing term", titlefont=12)

heatmap(x, y, vort_steady, xlabel="x", ylabel="y", title="Steady state vorticity", titlefont=12)

nothing

# The discrete Fourier transform of this vector fields are given as

g_steady_hat = rfft(g_steady)

vort_steady_hat = rfft(vort_steady)

nothing

# ## Testing the steady state

# We are now ready to test the steadyness of this vorticity field:

@testset "Check single-mode stable steady state" begin
    u_steady_hat = Hu_hat .* vort_steady_hat
    v_steady_hat = Hv_hat .* vort_steady_hat

    u_steady = irfft(u_steady_hat, N)
    v_steady = irfft(v_steady_hat, N)

    vort_steady_back = irfft(vort_steady_hat, N)

    wu_steady_hat = rfft(vort_steady_back .* u_steady)
    wv_steady_hat = rfft(vort_steady_back .* v_steady)

    rhs_steady = g_steady_hat .+ Delta_hat .* vort_steady_hat .- Dx_hat .* wu_steady_hat .- Dy_hat .* wv_steady_hat

    vort_steady_sol_hat = - g_steady_hat ./ Delta_hat
    vort_steady_sol_hat[1, 1] = 0.0im
    ## Vanishing bilinear term on one-mode steady state
    @test maximum(abs.(Dx_hat .* wu_steady_hat .+ Dy_hat .* wv_steady_hat)) ≤ √eps()
    ## Vanishing RHS on steady state
    @test maximum(abs.(rhs_steady)) ≤ √eps()
    ## Vanishing linear Stokes on steady state
    @test maximum(abs.(g_steady_hat .+ Delta_hat .* vort_steady_hat)) ≤ √eps()

    ## Steady state equation
    @test g_steady_hat ≈ - Delta_hat .* vort_steady_hat ≈ - Delta_hat .* vort_steady_hat .- Dx_hat .* wu_steady_hat .- Dy_hat .* wv_steady_hat

    ## Steady state solution
    @test vort_steady_sol_hat ≈ vort_steady_hat
end
nothing

# All seems good.


text:
using LevelSetAdvection
using Test

# test for Compliance
using Plots
gr()

# auxiliars
mean(x) = sum(x)/length(x)
mkpath("out")

# * == 1. Set the domain with Gridap
domain = (0,1,0,0.5)
partition = (180,90)
model = CartesianDiscreteModel(domain,partition)
topo  = get_grid_topology(model)
Ω     = Triangulation(model)

# and get some info
n      = num_cells(model) # num. elements
nv     = num_vertices(model)
dim    = num_dims(model) # 2

# extracts coordinates from the cell centers
xc = get_cellcenter_coordinates(Ω)
d_max = maximum(get_cell_diameter(Ω))

# * == 2. Set boundary conditions
labels = get_face_labeling(model)
add_tag_from_tags!(labels,"dirichlet",[1,3,7])
add_tag_from_tags!(labels,"load",8)

uD(x) = VectorValue(0.0,0.0)
g(x) = VectorValue(0.0,-10*(x[2]>=0.22)*(x[2]<=0.28))

order = 1
reffe = ReferenceFE(lagrangian,VectorValue{2,Float64},order)
V0 = TestFESpace(model,reffe, conformity=:H1, dirichlet_tags=["dirichlet"]) #  dirichlet_masks=[(true,false), (true,true)])

U  = TrialFESpace(V0,[uD])
degree = 2*order
Ω  = Triangulation(model)
dΩ = Measure(Ω,degree)
Γ  = BoundaryTriangulation(model,tags="load")
dΓ = Measure(Γ,degree)

# * == 3. Levelset initialization
#ϕ_(x) = min.((x[1] .- 0.25).^2 .+ (x[2] .- 0.25).^2 - 0.1^2,(x[1] .- 0.1).^2 .+ (x[2] .- 0.1).^2 - 0.08^2)
ϕ_(x) = -Signum.( sin.(4*2π*x[1]).*cos.(4*4π*x[2]) ,β=2)
ϕ = lazy_map(ϕ_,xc)

each_save  = 10
each_reinit = 2
max_iter   = 100

phi = vec(collect(get_array(ϕ))) ##! AVOID! Too slow
phi = ReinitHJ2d_update(topo,xc,phi,20,scheme="RK2",Δt=0.1*d_max)

# * == 4. Material setting
E = 10e9 * (phi.<=0) + 1e2*(phi.>0)
ν = 0.3

function E_to_C(E)

  λ = (E * ν)/((1+ν)*(1-2*ν))
  μ = E/(2*(1+ν))

  C1111 = 2.0*μ + λ
  C1122 = λ
  C1112 = 0.0
  C2222 = 2.0*μ + λ
  C2212 = 0.0
  C1212 = μ
  SymFourthOrderTensorValue(C1111,C1112,C1122,C1112,C1212,C2212,C1122,C2212,C2222)
end
C = lazy_map(E_to_C,E)

# Linear elasticity evaluation: First time
σ(u) = C⊙ε(u)
a(u,v) = ∫( ε(v) ⊙ σ(u) )dΩ
l(v) = ∫(v⋅g)dΓ

op = AffineFEOperator(a,l,U,V0)
uh = solve(op)

# * == 5. Compliance speed
compliance = [] # init

η = 1.0 # 0.2 # volume penalty
area = collect(get_array(∫(1.0)dΩ))
V(u) = σ(u) ⊙ ε(u)
Vc = collect(get_array(∫(V(uh))dΩ)) ./ area #! SLOW
# limiting speed near the boundary
#restric = @. exp(-abs(Vc)/(2*d_max))
#restric = @. restric*(restric>0.45)
#Vc .*= restric
Vc /= mean(Vc) #maximum(Vc)
println("[0] min,max(V) = ",minimum(Vc)," , ",maximum(Vc))
Δt = 0.2*d_max #/maximum(Vc) #<<<< does not allow 'maximum' over lazy array

writevtk(Ω,"out/elasticity_000",cellfields=["uh"=>uh,"epsi"=>ε(uh),"sigma"=>σ(uh)],celldata=["speed"=>Vc,"E"=>E,"phi"=>phi])

# * == 6. Optimization Loop

push!(compliance,sum(l(uh)))

#= only testing:
struct my_update <: Map end
evaluate!(cache,::my_update,x,y)=x-y
Update = my_update() =#

for k in 1:max_iter
  #global ϕ # for lazy_map
  global phi #? this is annoying
  ∇ϕ = upwind2d_step(topo,xc,phi,Vc.-η)

  #global phi = phi .- Δt*∇ϕ # Scope error in Julia
  #=
  for i in eachindex(phi)
    phi[i] -= Δt*∇ϕ[i] # this works
  end =#
  #ϕ = lazy_map(=,ϕ .- Δt.*∇ϕ) # not working
  #ϕ = lazy_map(-,Δt.*∇ϕ) #! ϕ is not updated! just local var!!
  phi = lazy_map(-,phi,Δt.*∇ϕ) #! ϕ not defined!!!! uses global ϕ

  # Reinitialization step
  if mod(k,each_reinit)==0
    phi = ReinitHJ2d_update(topo,xc,phi,10,scheme="RK2",Δt=0.1*d_max)
    #! HERE throws an error:
    # setindex! not defined for Gridap.Arrays.LazyArray (in LevelSetAdvection.jl:302)
  end

  global E = 10.0e9 * (phi.<=0) + 1e2*(phi.>0)
  global C = lazy_map(E_to_C,E)

  global op = AffineFEOperator(a,l,U,V0) # needed to update C into the problem
  global uh = solve(op)

  global Vc = collect(get_array(∫(V(uh))dΩ)) ./ area
  push!(compliance,sum(l(uh)))

  Vc /= mean(Vc) #maximum(Vc)
  println("[$(k)] min,max(V) = ",minimum(Vc)," , ",maximum(Vc))

  if mod(k,each_save)==0
    println("iter: ",k)
    writevtk(Ω,"out/elasticity_"*lpad(k,3,"0"),cellfields=["uh"=>uh,"epsi"=>ε(uh),"sigma"=>σ(uh)],celldata=["speed"=>Vc,"E"=>E,"phi"=>phi])
  end

  pp = plot(compliance, yaxis=:log10)
  display(pp)
  sleep(0.1)

end

max_stars_repo_name: joaquinmura/lsAdvect

max_issues_repo_name: joaquinmura/lsAdvect

max_forks_repo_name: joaquinmura/lsAdvect


text:
# Automatically generated using Clang.jl

const IVIDMM_MAJOR_VERSION = Int32(4)
const IVIDMM_MINOR_VERSION = Int32(0)
const IVIDMM_CLASS_SPEC_MAJOR_VERSION = Int32(3)
const IVIDMM_CLASS_SPEC_MINOR_VERSION = Int32(0)
const IVIDMM_DRIVER_VENDOR = "National Instruments"
const IVIDMM_DRIVER_DESCRIPTION = "IviDmm Class Driver"
const IVIDMM_VAL_MAX_TIME_INFINITE = IVI_VAL_MAX_TIME_INFINITE
const IVIDMM_VAL_MAX_TIME_IMMEDIATE = IVI_VAL_MAX_TIME_IMMEDIATE
const IVIDMM_ATTR_CACHE = IVI_ATTR_CACHE
const IVIDMM_ATTR_RANGE_CHECK = IVI_ATTR_RANGE_CHECK
const IVIDMM_ATTR_QUERY_INSTRUMENT_STATUS = IVI_ATTR_QUERY_INSTRUMENT_STATUS
const IVIDMM_ATTR_RECORD_COERCIONS = IVI_ATTR_RECORD_COERCIONS
const IVIDMM_ATTR_SIMULATE = IVI_ATTR_SIMULATE
const IVIDMM_ATTR_INTERCHANGE_CHECK = IVI_ATTR_INTERCHANGE_CHECK
const IVIDMM_ATTR_SPY = IVI_ATTR_SPY
const IVIDMM_ATTR_USE_SPECIFIC_SIMULATION = IVI_ATTR_USE_SPECIFIC_SIMULATION
const IVIDMM_ATTR_GROUP_CAPABILITIES = IVI_ATTR_GROUP_CAPABILITIES
const IVIDMM_ATTR_FUNCTION_CAPABILITIES = IVI_ATTR_FUNCTION_CAPABILITIES
const IVIDMM_ATTR_CLASS_DRIVER_PREFIX = IVI_ATTR_CLASS_DRIVER_PREFIX
const IVIDMM_ATTR_CLASS_DRIVER_VENDOR = IVI_ATTR_CLASS_DRIVER_VENDOR
const IVIDMM_ATTR_CLASS_DRIVER_DESCRIPTION = IVI_ATTR_CLASS_DRIVER_DESCRIPTION
const IVIDMM_ATTR_CLASS_DRIVER_CLASS_SPEC_MAJOR_VERSION = IVI_ATTR_CLASS_DRIVER_CLASS_SPEC_MAJOR_VERSION
const IVIDMM_ATTR_CLASS_DRIVER_CLASS_SPEC_MINOR_VERSION = IVI_ATTR_CLASS_DRIVER_CLASS_SPEC_MINOR_VERSION
const IVIDMM_ATTR_SPECIFIC_DRIVER_PREFIX = IVI_ATTR_SPECIFIC_DRIVER_PREFIX
const IVIDMM_ATTR_SPECIFIC_DRIVER_LOCATOR = IVI_ATTR_SPECIFIC_DRIVER_LOCATOR
const IVIDMM_ATTR_IO_RESOURCE_DESCRIPTOR = IVI_ATTR_IO_RESOURCE_DESCRIPTOR
const IVIDMM_ATTR_LOGICAL_NAME = IVI_ATTR_LOGICAL_NAME
const IVIDMM_ATTR_SPECIFIC_DRIVER_VENDOR = IVI_ATTR_SPECIFIC_DRIVER_VENDOR
const IVIDMM_ATTR_SPECIFIC_DRIVER_DESCRIPTION = IVI_ATTR_SPECIFIC_DRIVER_DESCRIPTION
const IVIDMM_ATTR_SPECIFIC_DRIVER_CLASS_SPEC_MAJOR_VERSION = IVI_ATTR_SPECIFIC_DRIVER_CLASS_SPEC_MAJOR_VERSION
const IVIDMM_ATTR_SPECIFIC_DRIVER_CLASS_SPEC_MINOR_VERSION = IVI_ATTR_SPECIFIC_DRIVER_CLASS_SPEC_MINOR_VERSION
const IVIDMM_ATTR_INSTRUMENT_FIRMWARE_REVISION = IVI_ATTR_INSTRUMENT_FIRMWARE_REVISION
const IVIDMM_ATTR_INSTRUMENT_MANUFACTURER = IVI_ATTR_INSTRUMENT_MANUFACTURER
const IVIDMM_ATTR_INSTRUMENT_MODEL = IVI_ATTR_INSTRUMENT_MODEL
const IVIDMM_ATTR_SUPPORTED_INSTRUMENT_MODELS = IVI_ATTR_SUPPORTED_INSTRUMENT_MODELS
const IVIDMM_ATTR_CLASS_DRIVER_REVISION = IVI_ATTR_CLASS_DRIVER_REVISION
const IVIDMM_ATTR_SPECIFIC_DRIVER_REVISION = IVI_ATTR_SPECIFIC_DRIVER_REVISION
const IVIDMM_ATTR_DRIVER_SETUP = IVI_ATTR_DRIVER_SETUP
const IVIDMM_ATTR_FUNCTION = IVI_CLASS_PUBLIC_ATTR_BASE + Int32(1)
const IVIDMM_ATTR_RANGE = IVI_CLASS_PUBLIC_ATTR_BASE + Int32(2)
const IVIDMM_ATTR_RESOLUTION_ABSOLUTE = IVI_CLASS_PUBLIC_ATTR_BASE + Int32(8)
const IVIDMM_ATTR_TRIGGER_SOURCE = IVI_CLASS_PUBLIC_ATTR_BASE + Int32(4)
const IVIDMM_ATTR_TRIGGER_DELAY = IVI_CLASS_PUBLIC_ATTR_BASE + Int32(5)
const IVIDMM_ATTR_AC_MIN_FREQ = IVI_CLASS_PUBLIC_ATTR_BASE + Int32(6)
const IVIDMM_ATTR_AC_MAX_FREQ = IVI_CLASS_PUBLIC_ATTR_BASE + Int32(7)
const IVIDMM_ATTR_FREQ_VOLTAGE_RANGE = IVI_CLASS_PUBLIC_ATTR_BASE + Int32(101)
const IVIDMM_ATTR_TEMP_TRANSDUCER_TYPE = IVI_CLASS_PUBLIC_ATTR_BASE + Int32(201)
const IVIDMM_ATTR_TEMP_TC_TYPE = IVI_CLASS_PUBLIC_ATTR_BASE + Int32(231)
const IVIDMM_ATTR_TEMP_TC_REF_JUNC_TYPE = IVI_CLASS_PUBLIC_ATTR_BASE + Int32(232)
const IVIDMM_ATTR_TEMP_TC_FIXED_REF_JUNC = IVI_CLASS_PUBLIC_ATTR_BASE + Int32(233)
const IVIDMM_ATTR_TEMP_RTD_ALPHA = IVI_CLASS_PUBLIC_ATTR_BASE + Int32(241)
const IVIDMM_ATTR_TEMP_RTD_RES = IVI_CLASS_PUBLIC_ATTR_BASE + Int32(242)
const IVIDMM_ATTR_TEMP_THERMISTOR_RES = IVI_CLASS_PUBLIC_ATTR_BASE + Int32(251)
const IVIDMM_ATTR_SAMPLE_COUNT = IVI_CLASS_PUBLIC_ATTR_BASE + Int32(301)
const IVIDMM_ATTR_SAMPLE_TRIGGER = IVI_CLASS_PUBLIC_ATTR_BASE + Int32(302)
const IVIDMM_ATTR_SAMPLE_INTERVAL = IVI_CLASS_PUBLIC_ATTR_BASE + Int32(303)
const IVIDMM_ATTR_TRIGGER_COUNT = IVI_CLASS_PUBLIC_ATTR_BASE + Int32(304)
const IVIDMM_ATTR_MEAS_COMPLETE_DEST = IVI_CLASS_PUBLIC_ATTR_BASE + Int32(305)
const IVIDMM_ATTR_TRIGGER_SLOPE = IVI_CLASS_PUBLIC_ATTR_BASE + Int32(334)
const IVIDMM_ATTR_APERTURE_TIME = IVI_CLASS_PUBLIC_ATTR_BASE + Int32(321)
const IVIDMM_ATTR_APERTURE_TIME_UNITS = IVI_CLASS_PUBLIC_ATTR_BASE + Int32(322)
const IVIDMM_ATTR_AUTO_RANGE_VALUE = IVI_CLASS_PUBLIC_ATTR_BASE + Int32(331)
const IVIDMM_ATTR_AUTO_ZERO = IVI_CLASS_PUBLIC_ATTR_BASE + Int32(332)
const IVIDMM_ATTR_POWERLINE_FREQ = IVI_CLASS_PUBLIC_ATTR_BASE + Int32(333)
const IVIDMM_VAL_DC_VOLTS = Int32(1)
const IVIDMM_VAL_AC_VOLTS = Int32(2)
const IVIDMM_VAL_DC_CURRENT = Int32(3)
const IVIDMM_VAL_AC_CURRENT = Int32(4)
const IVIDMM_VAL_2_WIRE_RES = Int32(5)
const IVIDMM_VAL_4_WIRE_RES = Int32(101)
const IVIDMM_VAL_AC_PLUS_DC_VOLTS = Int32(106)
const IVIDMM_VAL_AC_PLUS_DC_CURRENT = Int32(107)
const IVIDMM_VAL_FREQ = Int32(104)
const IVIDMM_VAL_PERIOD = Int32(105)
const IVIDMM_VAL_TEMPERATURE = Int32(108)
const IVIDMM_VAL_FUNC_CLASS_EXT_BASE = Int32(500)
const IVIDMM_VAL_FUNC_SPECIFIC_EXT_BASE = Int32(1000)
const IVIDMM_VAL_AUTO_RANGE_ON = -1.0
const IVIDMM_VAL_AUTO_RANGE_OFF = -2.0
const IVIDMM_VAL_AUTO_RANGE_ONCE = -3.0
const IVIDMM_VAL_RANGE_CLASS_EXT_BASE = -100.0
const IVIDMM_VAL_RANGE_SPECIFIC_EXT_BASE = -1000.0
const IVIDMM_VAL_FREQ_VOLT_RANGE_CLASS_EXT_BASE = -(Float64(100.0))
const IVIDMM_VAL_FREQ_VOLT_RANGE_SPECIFIC_EXT_BASE = -(Float64(1000.0))
const IVIDMM_VAL_AUTO_DELAY_ON = -(Float64(1.0))
const IVIDMM_VAL_AUTO_DELAY_OFF = -(Float64(2.0))
const IVIDMM_VAL_TRIGGER_DELAY_CLASS_EXT_BASE = -(Float64(100.0))
const IVIDMM_VAL_TRIGGER_DELAY_SPECIFIC_EXT_BASE = -(Float64(1000.0))
const IVIDMM_VAL_IMMEDIATE = Int32(1)
const IVIDMM_VAL_EXTERNAL = Int32(2)
const IVIDMM_VAL_SOFTWARE_TRIG = Int32(3)
const IVIDMM_VAL_TTL0 = Int32(111)
const IVIDMM_VAL_TTL1 = Int32(112)
const IVIDMM_VAL_TTL2 = Int32(113)
const IVIDMM_VAL_TTL3 = Int32(114)
const IVIDMM_VAL_TTL4 = Int32(115)
const IVIDMM_VAL_TTL5 = Int32(116)
const IVIDMM_VAL_TTL6 = Int32(117)
const IVIDMM_VAL_TTL7 = Int32(118)
const IVIDMM_VAL_ECL0 = Int32(119)
const IVIDMM_VAL_ECL1 = Int32(120)
const IVIDMM_VAL_PXI_STAR = Int32(131)
const IVIDMM_VAL_RTSI_0 = Int32(140)
const IVIDMM_VAL_RTSI_1 = Int32(141)
const IVIDMM_VAL_RTSI_2 = Int32(142)
const IVIDMM_VAL_RTSI_3 = Int32(143)
const IVIDMM_VAL_RTSI_4 = Int32(144)
const IVIDMM_VAL_RTSI_5 = Int32(145)
const IVIDMM_VAL_RTSI_6 = Int32(146)
const IVIDMM_VAL_TRIGGER_SOURCE_CLASS_EXT_BASE = Int32(500)
const IVIDMM_VAL_TRIGGER_SOURCE_SPECIFIC_EXT_BASE = Int32(1000)
const IVIDMM_VAL_THERMOCOUPLE = Int32(1)
const IVIDMM_VAL_THERMISTOR = Int32(2)
const IVIDMM_VAL_2_WIRE_RTD = Int32(3)
const IVIDMM_VAL_4_WIRE_RTD = Int32(4)
const IVIDMM_VAL_TRANSDUCER_CLASS_EXT_BASE = Int32(100)
const IVIDMM_VAL_TRANSDUCER_SPECIFIC_EXT_BASE = Int32(1000)
const IVIDMM_VAL_TEMP_REF_JUNC_INTERNAL = Int32(1)
const IVIDMM_VAL_TEMP_REF_JUNC_FIXED = Int32(2)
const IVIDMM_VAL_TEMP_REF_JUNC_CLASS_EXT_BASE = Int32(100)
const IVIDMM_VAL_TEMP_REF_JUNC_SPECIFIC_EXT_BASE = Int32(1000)
const IVIDMM_VAL_TEMP_TC_B = Int32(1)
const IVIDMM_VAL_TEMP_TC_C = Int32(2)
const IVIDMM_VAL_TEMP_TC_D = Int32(3)
const IVIDMM_VAL_TEMP_TC_E = Int32(4)
const IVIDMM_VAL_TEMP_TC_G = Int32(5)
const IVIDMM_VAL_TEMP_TC_J = Int32(6)
const IVIDMM_VAL_TEMP_TC_K = Int32(7)
const IVIDMM_VAL_TEMP_TC_N = Int32(8)
const IVIDMM_VAL_TEMP_TC_R = Int32(9)
const IVIDMM_VAL_TEMP_TC_S = Int32(10)
const IVIDMM_VAL_TEMP_TC_T = Int32(11)
const IVIDMM_VAL_TEMP_TC_U = Int32(12)
const IVIDMM_VAL_TEMP_TC_V = Int32(13)
const IVIDMM_VAL_TEMP_TC_TYPE_CLASS_EXT_BASE = Int32(100)
const IVIDMM_VAL_TEMP_TC_TYPE_SPECIFIC_EXT_BASE = Int32(1000)
const IVIDMM_VAL_NONE = -(Int32(1))
const IVIDMM_VAL_INTERVAL = Int32(10)
const IVIDMM_VAL_POSITIVE = Int32(0)
const IVIDMM_VAL_NEGATIVE = Int32(1)
const IVIDMM_VAL_TRIGGER_SLOPE_CLASS_EXT_BASE = Int32(100)
const IVIDMM_VAL_TRIGGER_SLOPE_SPECIFIC_EXT_BASE = Int32(1000)
const IVIDMM_VAL_SECONDS = Int32(0)
const IVIDMM_VAL_POWER_LINE_CYCLES = Int32(1)
const IVIDMM_VAL_AUTO_ZERO_OFF = Int32(0)
const IVIDMM_VAL_AUTO_ZERO_ON = Int32(1)
const IVIDMM_VAL_AUTO_ZERO_ONCE = Int32(2)
const IVIDMM_VAL_AUTO_ZERO_CLASS_EXT_BASE = Int32(100)
const IVIDMM_VAL_AUTO_ZERO_SPECIFIC_EXT_BASE = Int32(1000)
const IVIDMM_WARN_OVER_RANGE = IVI_CLASS_WARN_BASE + 1
const IVIDMM_ERROR_TRIGGER_NOT_SOFTWARE = IVI_CROSS_CLASS_ERROR_BASE + 1
const IVIDMM_ERROR_MAX_TIME_EXCEEDED = IVI_CLASS_ERROR_BASE + 3
const IVIDMM_WRNMSG_OVER_RANGE = "Over range."
const IVIDMM_ERRMSG_MAX_TIME_EXCEEDED = "Max Time exceeded before "
const IVIDMM_ERRMSG_TRIGGER_NOT_SOFTWARE = "The trigger source is not set to "

# Skipping MacroDefinition: IVIDMM_ERROR_CODES_AND_MSGS { IVIDMM_WARN_OVER_RANGE , IVIDMM_WRNMSG_OVER_RANGE } , { IVIDMM_ERROR_MAX_TIME_EXCEEDED , IVIDMM_ERRMSG_MAX_TIME_EXCEEDED } , { IVIDMM_ERROR_TRIGGER_NOT_SOFTWARE , IVIDMM_ERRMSG_TRIGGER_NOT_SOFTWARE }


text:
using Pkg
using CodecZlib
import RegistryTools

function gitcmd(config, path)
    return RegistryTools.gitcmd(path, config.gitconfig)
end

function get_local_registry_dir(config)
    dir = joinpath(config.git_clones_dir, "registry")
    return dir
end

function get_local_package_dir(config, uuid)
    dir = joinpath(config.git_clones_dir, "packages", uuid)
    return dir
end

function get_registries(config, server::GitStorageServer)
    repo = server.url
    registry_dir = get_local_registry_dir(config)
    git = gitcmd(config, registry_dir)
    if !isdir(registry_dir)
        mkpath(registry_dir)
        run(`$git clone $(repo) .`)
    else
        run(`$git pull --quiet`)
    end
    registry = Pkg.TOML.parsefile(joinpath(registry_dir, "Registry.toml"))
    server.uuid = registry["uuid"]
    hash = readchomp(`$git rev-parse --verify HEAD:`)
    return Dict(registry["uuid"] => hash)
end

function get_resource_from_storage_server!(config, server::GitStorageServer,
                                           resource, io)
    parts = split(resource, "/", keepempty = false)
    registry_dir = get_local_registry_dir(config)
    if parts[1] == "registry"
        git = gitcmd(config, registry_dir)
        uuid = parts[2]
        hash = parts[3]
        uuid != server.uuid && return false
    elseif parts[1] == "package"
        uuid = parts[2]
        hash = parts[3]
        # TODO: Cache Registry.toml in memory.
        registry = Pkg.TOML.parsefile(joinpath(registry_dir, "Registry.toml"))
        haskey(registry["packages"], uuid) || return false
        package_dir = joinpath(registry_dir, registry["packages"][uuid]["path"])
        package = Pkg.TOML.parsefile(joinpath(package_dir, "Package.toml"))
        repo = package["repo"]
        local_package_dir = get_local_package_dir(config, uuid)
        git = gitcmd(config, local_package_dir)
        if !isdir(local_package_dir)
            @info "Cloning package" repo=repo Dates.now()
            mkpath(local_package_dir)
            run(`$git clone --mirror $(repo) .`)
        else
            # If hash is not available, update git repo.
            try
                run(`$git rev-parse --verify --quiet $(hash)^\{tree\}`)
            catch
                @info "Hash not available, updating from remote" repo=repo hash=hash Dates.now()
                run(`$git remote update`)
            end
            # Still not there? Nothing to do about it.
            try
                run(`$git rev-parse --verify --quiet $(hash)^\{tree\}`)
            catch
                @error "Hash still not available after remote update" repo=repo hash=hash Dates.now()
                close(io)
                return false
            end
        end
    elseif parts[1] == "artifact"
        @info "No support for artifacts" Dates.now()
        return false
    else
        @info "Unknown resource $(parts[1])" Dates.now()
        return false
    end

    # Do not allow git on windows to convert line endings in the tarball.
    tar = read(`$git -c core.autocrlf=false archive $(hash)`)
    gzip = transcode(GzipCompressor, tar)
    write(io, gzip)
    close(io)
    return true
end


text:
import Base.show, Base.==, Base.adjoint, Base.*, Base.eltype
import Base.getindex
import LightXML.name

export SimpleGraph, IntGraph, StringGraph
export show, NV, NE, has, typ, fastN!, name, get_edge
export vlist, elist, neighbors, getindex, deg, deg_hist

"""
The `SimpleGraph` type represents a simple graph; that is, an
undirected graph with no loops and no multiple edges.

Use `SimpleGraph()` to create a new graph in which the vertices may be
`Any` type. Use `SimpleGraph{T}()` to create a new graph in which the
vertices are of type `T`. See `IntGraph` and `StringGraph` as special
cases.
"""
mutable struct SimpleGraph{T} <: AbstractSimpleGraph
    V::Set{T}          # Vertex set
    E::Set{Tuple{T,T}} # Edge set
    N::Dict{T,Set{T}}  # Optional neighbor sets
    Nflag::Bool        # Tells if N is used or not (default on)
    cache::Dict{Symbol,Any}   # save previous expensive results
    cache_flag::Bool          # decide if we use cache or no
    function SimpleGraph{T}(Nflag::Bool = true) where {T}
        V = Set{T}()
        E = Set{Tuple{T,T}}()
        N = Dict{T,Set{T}}()
        cache = Dict{Symbol,Any}()
        G = new(V, E, N, Nflag, cache, true)
    end
end

"""
`name(G)` returns the graph's name.

`name(G,str)` assigns `str` to be the graph's name. If `str` is
empty, then the name is set to the default `SimpleGraph{T}` where
`T` is the vertex type.
"""
function name(G::SimpleGraph)
    if cache_check(G, :name)
        return cache_recall(G, :name)
    end
    return "SimpleGraph{$(eltype(G))}"
end

function name(G::SimpleGraph, the_name::String)
    G.cache[:name] = the_name
    if length(the_name) == 0
        cache_clear(G, :name)
    end
    nothing
end

function show(io::IO, G::SimpleGraph)
    suffix = " (n=$(NV(G)), m=$(NE(G)))"
    print(io, name(G) * suffix)
end

# Default constructor uses Any type vertices
SimpleGraph(Nflag::Bool = true) = SimpleGraph{Any}(Nflag)

# A StringGraph has vertices of type String.
"""
A `StringGraph` is a `SimpleGraph` whose vertices are of type
`String`.

When constructed with `StringGraph()` creates an empty
`SimpleGraph{String}`.

When invoked as `StringGraph(file::AbstractString)` opens the named
file for parsing as a graph. Each line of the file should contain one
or two tokens separated by white space. If the line contains a single
token, we add that token as a vertex. If the line contains two (or
more) tokens, then the first two tokens are taken as vertex names and
(assuming they are different) the corresponding edge is created. Any
extra tokens on the line are ignored. Lines that begin with a # are
ignored.
"""
StringGraph() = SimpleGraph{String}()

function StringGraph(file::AbstractString)
    G = StringGraph()
    load!(G, file)
    return G
end

# Helper function for StrinGraph(file), and can be used to add
# vertices and edges to a graph (assuming its vertex type can
# accomodate strings).
function load!(G::SimpleGraph, file::AbstractString)
    f = open(file, "r")
    while (~eof(f))
        line = chomp(readline(f))
        tokens = split(line)

        if (length(tokens) == 0)
            continue
        end

        if (tokens[1][1] == '#')
            continue
        end

        add!(G, tokens[1])
        if length(tokens) > 1
            add!(G, tokens[1], tokens[2])
        end
    end
end

"""
`IntGraph()` creates a new `SimpleGraph` whose vertices are of type
`Int`. Called as `IntGraph(n::Int)` prepopulates the vertex set with
vertices `1:n`.

`IntGraph(A)` where `A` is an adjacency matrix creates a graph for which
`A` is the adjacency matrix.
"""
IntGraph() = SimpleGraph{Int}()

# With a postive integer argument, adds 1:n as vertex set, but no
# edges.
function IntGraph(n::Int)
    G = IntGraph()
    for v = 1:n
        add!(G, v)
    end
    return G
end

function IntGraph(A::AbstractMatrix)
    r, c = size(A)
    @assert r == c "Matrix must be square"
    @assert A == A' "Matrix must be symmetric"
    G = IntGraph(r)
    for i = 1:r-1
        for j = i+1:r
            if A[i, j] != 0
                add!(G, i, j)
            end
        end
    end
    return G
end

"""
`SimpleGraph(A)` where `A` is a matrix creates a graph with vertex set
`1:n` where `A` is an `n`-by-`n` symmetric matrix specifying the graph's
adjacency matrix.
"""
SimpleGraph(A::AbstractMatrix) = IntGraph(A)

"""
`eltype(G::SimpleGraph)`, `eltype(G::SimpleDigraph)`, `eltype(G::SimpleHypergraph)`

Returns the data type of the vertices this graph may hold.
For example, if `G=IntGraph()` then this returns `Int64`.`
"""
eltype(G::SimpleGraph{T}) where {T} = T

@deprecate vertex_type eltype

# number of vertices and edges

"""
`NV(G)` returns the number of vertices in `G`.
"""
NV(G::AbstractSimpleGraph) = length(G.V)

"""
`NE(G)` returns the number of edges in `G`.
"""
NE(G::SimpleGraph) = length(G.E)

"""
`has(G,v)` returns `true` iff `v` is a vertex of `G`.

`has(G,v,w)` returns `true` iff `(v,w)` is an edge of `G`.
"""
has(G::AbstractSimpleGraph, v) = in(v, G.V)
has(G::SimpleGraph, v, w) = in((v, w), G.E) || in((w, v), G.E)

"""
`get_edge(G,u,v)` returns either `(u,v)` or `(v,u)` matching  how
the edge joining `u` and `v` is stored in the edge set of `G`.
An error is thrown if `u` and `v` are not adjacent vertices of `G`.
"""
function get_edge(G::SimpleGraph{T}, u, v)::Tuple{T,T} where {T}
    if !has(G, u, v)
        error("($u,$v) is not an edge of this graph")
    end
    if in((u, v), G.E)
        return u, v
    end
    return v, u
end

# fastN(G,true) creates an additional data structure to speed up
# neighborhood lookups.

"""
`fastN!(G,flg=true)` is used to turn on (or off) fast neighborhood
lookup in graphs. Switching this off decreases the size of the data
structure holding the graph, but slows down look up of edges.

**Note**: Fast neighborhood look up is on by default.
"""
function fastN!(G::SimpleGraph{T}, flg::Bool = true) where {T}
    # if no change, do nothing
    if flg == G.Nflag
        return flg
    end

    # if setting the flag to false, erase the G.N data structure;
    # otherwise build it.
    if flg == false
        G.N = Dict{T,Set{T}}()
        sizehint!(G.N, NV(G))
    else
        # build the G.N structure.
        # start with empty sets for each vertex
        for v in G.V
            G.N[v] = Set{T}()
        end
        # now iterate over the edge set and populate G.N sets
        for e in G.E
            v, w = e
            push!(G.N[v], w)
            push!(G.N[w], v)
        end
    end
    G.Nflag = flg
    return flg
end

# Create a mapping between G.V and 1:n. This is not exposed outside
# this module; it's a helper function used by other functions. This
# has been crafted to work with either SimpleGraph or SimpleDigraph
# arguments.
function vertex2idx(G::AbstractSimpleGraph)
    T = eltype(G)
    d = Dict{T,Int}()
    V = vlist(G)
    n = NV(G)

    for k = 1:n
        d[V[k]] = k
    end

    return d
end

# get the vertices as a (sorted if possible) list
"""
`vlist(G)` returns the vertices of `G` as a list (array).
"""
function vlist(G::AbstractSimpleGraph)
    result = collect(G.V)
    try
        sort!(result)
    catch
    end
    return result
end

"""
`elist(G)` returns the edges of `G` as a list (array).
"""
function elist(G::SimpleGraph)
    result = collect(G.E)
    try
        sort!(result)
    catch
    end
    return result
end

# Get the neighbors of a vertex
"""
`neighbors(G,v)` returns a list of the neighbors of `v`.

May also be invoked as `G[v]`.
"""
function neighbors(G::SimpleGraph{T}, v) where {T}
    if ~has(G, v)
        error("Graph does not contain requested vertex")
    end

    if G.Nflag
        N = collect(G.N[v])

    else
        N = T[]
        for w in G.V
            if has(G, v, w)
                append!(N, [w])
            end
        end
    end
    try
        sort!(N)
    catch
    end
    return N
end

# Here is another way to access the neighbors of a vertex: G[v]
getindex(G::SimpleGraph, v) = neighbors(G, v)

# And here's a getindex way to check for edges: G[u,v] is a shortcut
# for has(G,u,v).
getindex(G::SimpleGraph, v, w) = has(G, v, w)

# Degree of a vertex
"""
`deg(G,v)` gives the degree of `v` in `G`.

`deg(G)` gives the degree sequence (sorted).
"""
function deg(G::SimpleGraph, v)
    if ~has(G, v)
        error("Graph does not contain requested vertex")
    end
    if G.Nflag
        return length(G.N[v])
    end
    return length(G[v])
end

# Degree sequence
function deg(G::SimpleGraph{T}) where {T}
    if G.Nflag
        ds = [deg(G, v) for v in G.V]
    else
        dd = Dict{T,Int}()
        for v in G.V
            dd[v] = 0
        end
        for e in G.E
            v, w = e
            dd[v] += 1
            dd[w] += 1
        end
        ds = collect(values(dd))
    end
    sort!(ds, lt = >)
    return ds
end

# Report how many vertices we have each possible degree.
# If G has n vertices, this returns an n-long vector whose
# k'th entry is the number of vertices of degree k-1.

"""
`deg_hist(G)` gives a tally of how many vertices of each degree are
present in the graph. Because Julia arrays are 1-based, the indexing
is a bit off. Specifically, entry `k` in the returned array is the
number of vertices of degree `k-1`.
"""
function deg_hist(G::SimpleGraph{T}) where {T}
    n = NV(G)
    degs = deg(G)
    result = zeros(Int, n)
    for d in degs
        result[d+1] += 1
    end
    return result
end

import Base.hash

function hash(G::SimpleGraph, h::UInt64 = UInt64(0))
    return hash(G.V, h) + hash(G.E, h)
end

simplify(G::SimpleGraph) = deepcopy(G)


text:
"""
    abstract type BoundaryTriangulation{Dc,Dp} <: Triangulation{Dc,Dp}
"""
abstract type BoundaryTriangulation{Dc,Dp} <: Triangulation{Dc,Dp} end

"""
    get_volume_triangulation(trian::BoundaryTriangulation)
"""
function get_volume_triangulation(trian::BoundaryTriangulation)
  @abstractmethod
end

"""
"""
function get_face_to_cell(trian::BoundaryTriangulation)
  @abstractmethod
end

"""
"""
function get_face_to_lface(trian::BoundaryTriangulation)
  @abstractmethod
end

"""
"""
function get_face_to_cell_map(trian::BoundaryTriangulation)
  @abstractmethod
end

"""
"""
function get_normal_vector(trian::BoundaryTriangulation)
  @abstractmethod
end

"""
"""
function get_face_to_face(trian::BoundaryTriangulation)
  @abstractmethod
end

"""
"""
function get_cell_around(trian::BoundaryTriangulation)
  @abstractmethod
end

# Tester

"""
    test_boundary_triangulation(trian::BoundaryTriangulation)
"""
function test_boundary_triangulation(trian::BoundaryTriangulation)
  test_triangulation(trian)
  @test isa(get_volume_triangulation(trian),Triangulation)
  @test isa(get_face_to_cell(trian),AbstractArray{<:Integer})
  @test isa(get_face_to_lface(trian),AbstractArray{<:Integer})
  @test isa(get_face_to_cell_map(trian),AbstractArray{<:Field})
  @test isa(get_normal_vector(trian),CellField)
end

# Default API

function restrict(f::AbstractArray, trian::BoundaryTriangulation)
  compose_field_arrays(reindex(f,trian), get_face_to_cell_map(trian))
end

function get_cell_id(trian::BoundaryTriangulation)
  get_face_to_cell(trian)
end

# Constructors

"""
    BoundaryTriangulation(model::DiscreteModel,face_to_mask::Vector{Bool})
    BoundaryTriangulation(model::DiscreteModel)
"""
function BoundaryTriangulation(model::DiscreteModel,face_to_mask::Vector{Bool},icell_around)
  GenericBoundaryTriangulation(model,face_to_mask,icell_around)
end

function BoundaryTriangulation(model::DiscreteModel,face_to_mask::Vector{Bool})
  icell_around = 1
  BoundaryTriangulation(model,face_to_mask,icell_around)
end

function BoundaryTriangulation(model::DiscreteModel)
  topo = get_grid_topology(model)
  D = num_cell_dims(model)
  face_to_mask = collect(Bool,get_isboundary_face(topo,D-1))
  BoundaryTriangulation(model,face_to_mask)
end

"""
    BoundaryTriangulation(model::DiscreteModel,tags::Vector{Int})
    BoundaryTriangulation(model::DiscreteModel,tags::Vector{String})
    BoundaryTriangulation(model::DiscreteModel,tag::Int)
    BoundaryTriangulation(model::DiscreteModel,tag::String)
"""
function BoundaryTriangulation(model::DiscreteModel,tags)
  labeling = get_face_labeling(model)
  BoundaryTriangulation(model,labeling,tags)
end

"""
    BoundaryTriangulation(model::DiscreteModel,labeling::FaceLabeling,tags::Vector{Int})
    BoundaryTriangulation(model::DiscreteModel,labeling::FaceLabeling,tags::Vector{String})
    BoundaryTriangulation(model::DiscreteModel,labeling::FaceLabeling,tag::Int)
    BoundaryTriangulation(model::DiscreteModel,labeling::FaceLabeling,tag::String)
"""
function BoundaryTriangulation(model::DiscreteModel,labeling::FaceLabeling,tags)
  D = num_cell_dims(model)
  face_to_mask = get_face_mask(labeling,tags,D-1)
  BoundaryTriangulation(model,face_to_mask)
end

#function BoundaryTriangulation(model::DiscreteModel,names::Vector{String})
#  labeling = get_face_labeling(model)
#  tags = get_tags_from_names(labeling,names)
#  BoundaryTriangulation(model,tags)
#end
#
#function BoundaryTriangulation(model::DiscreteModel,tag::Union{Int,String})
#  tags = [tag,]
#  BoundaryTriangulation(model,tags)
#end
#
#function _convert_to_face_to_masks(labeling,tags)
#  D = num_cell_dims(model)
#  face_to_mask = get_face_mask(labeling,tags,D-1)
#end


text:
""" Supports conversion from PetriNets to Catalyst ReactionSystems

This provides access to the parameter estimation, optimization, and sensitivity
tooling provided in the Catalyst library
"""

module CatalystInterop
  using AlgebraicPetri
  using Catlab.CategoricalAlgebra
  using ...Catalyst
  using ...Catalyst.Symbolics: scalarize
  import ...Catalyst: ReactionSystem

  counter(a) = [count(==(i),a) for i in unique(a)]

  """ Convert a general PetriNet to a ReactionSystem

  This conversion forgets any labels or rates provided, and converts all
  parameters and variables into symbols. It does preserve the ordering of
  transitions and states though (Transition 1 has a rate of k[1], state 1 has a
  concentration of S[1])
  """
  function ReactionSystem(pn::AbstractPetriNet)
    @parameters t k[1:nt(pn)]
    @variables S[collect(1:ns(pn))](t)

    rxs = map(1:nt(pn)) do t
      inpts = pn[incident(pn, t, :it),:is]
      otpts = pn[incident(pn, t, :ot),:os]
      in_count = collect(counter(inpts))
      ot_count = collect(counter(otpts))
      Reaction(k[t], [S[i] for i in unique(inpts)],
                     [S[o] for o in unique(otpts)],
                     in_count, ot_count)
    end

    ReactionSystem(rxs, t, scalarize(S), scalarize(k))
  end
end

max_stars_repo_name: AlgebraicJulia/AlgebraicPetri

max_issues_repo_name: AlgebraicJulia/AlgebraicPetri

max_forks_repo_name: AlgebraicJulia/AlgebraicPetri


text:
module SpawnFacing

using ..Ahorn, Maple

const placements = Ahorn.PlacementDict(
    "Spawn Facing" => Ahorn.EntityPlacement(
        Maple.SpawnFacingTrigger,
        "rectangle"
    )
)

function Ahorn.editingOptions(trigger::Maple.SpawnFacingTrigger)
    return Dict{String, Any}(
        "facing" => sort(Maple.spawn_facing_trigger_facings)
    )
end

end