text:
using Test
using Polynomials, DataFrames

push!(LOAD_PATH, joinpath(pwd(), "src")) # run from main dir
using Controlz

const s = TransferFunction([1, 0], [1])

@testset "testing TransferFunctions" begin
    ###
    # constructors
    ###
    @test TransferFunction([1], [3, 2], 0.0) == TransferFunction([1], [3, 2])
    tf = TransferFunction([1], [3, 2], 3.2)
    @test tf.time_delay == 3.2
    @test tf.numerator == Polynomial([1], :s)
    @test tf.denominator == Polynomial([2, 3], :s) # reverse order
    # special constructors
    @test first_order_system(5.0, 8.0) == 5 / (8 * s + 1)
    @test second_order_system(5.0, 8.0, 0.1) == 5 / (8^2 * s^2 + 2*8*0.1*s + 1)

    ###
    #  isapprox
    ###
    g1 = 4 * s / (2 * s + 1)
    g2 = 8 * s / (4 * s + 2)
    @test isapprox(g1, g2)
    g2 = 8 * s / (4 * s + 2) * exp(-0.01*s)
    @test ! isapprox(g1, g2)
    @test isapprox((s^2 + 1) / (3 * s + 3), 1/3*(s^2+1)/(s+1))
    @test ! isapprox(s^2 + 1 / (3 * s + 3), s*1/3*(s^2+1)/(s+1))
    @test isapprox((5 * s + 1) / (s ^ 2 + 4 * s + 5), TransferFunction([5, 1], [1, 4, 5]))
    @test isapprox((5 * s + 1) / (s ^ 2 + 4 * s + 5), zeros_poles_k([-1/5], [-2 + im, -2 - im], 5.0, time_delay=0.0))

    ###
    # multiply
    ###
    g1 = TransferFunction([1], [3, 1])
    g2 = TransferFunction([2], [5, 1])
    @test isapprox(zero_frequency_gain(g1 * g2), 2.0)
    @test isapprox(g1 * g2, 2/(3*s+1)/(5*s+1))
    @test isapprox(g1 * g2, zeros_poles_k([], [-1/3, -1/5], 2/15))
    g1 = TransferFunction([2], [4, 1], 4.0)
    g2 = TransferFunction([1], [2, 1], 2.0)
    @test isapprox(g1 * g2, TransferFunction([2], [8, 6, 1], 6.0))
    g = TransferFunction([2], [3, 4])
    @test isapprox(5 * g, TransferFunction([10], [3, 4]))

    # gain
    @test isapprox(zero_frequency_gain(5 / (s + 1)), 5.0)
    @test isapprox(zero_frequency_gain(5 / (s + 1) / (4 * s + 2)), 5.0 / 2)
    @test isapprox(zero_frequency_gain(5 * s / (s + 1) / s), 5.0)

    ###
    #  add
    ###
    g1 = 3 / (s + 2)
    g2 = 1 / (s + 4)
    @test isapprox(g1 + g2, (4*s+14) / (s^2 + 6*s + 8))

    ###
    #   subtract
    ###
    @test isapprox(-s, -1.0 * s)
    g = 3 / (s+2) - 2 / (s+2)
    @test isapprox(pole_zero_cancellation(g), 1/(s+2))
    @test isapprox(g, 1/(s+2))

    ###
    #  divide
    ###
    g1 = TransferFunction([2], [4, 1], 4.0)
    g2 = TransferFunction([1], [2, 1], 2.0)
    @test isapprox(g1 / g2, TransferFunction([4, 2], [4, 1], 2.0))
    @test isapprox(g1 / 2, TransferFunction([1], [4, 1], 4.0))
    @test isapprox(5 / (s+3), TransferFunction([5], [1, 3]))
    @test isapprox(33 / (2*s+3), TransferFunction([33], [2, 3]))
    @test isapprox((s+3) / ((s - 2) * (2*s-3)), TransferFunction([1, 3], [2, -7, 6]))

    ###
    #  power
    ###
    @test isapprox(s ^ 2, s * s)
    @test isapprox(s ^ 4, s * s * s * s)
    @test isapprox(s ^ 1, s)
    g = TransferFunction([2], [3, 4])
    @test isapprox(g * s ^ 0, g)
    @test_throws MethodError s ^ -1

    ###
    #  tf algebra
    ###
    τ = 1.4
    K = 2.8
    Kc = 3.3
    Kd = 4.5
    τI = 9.0
    τd = 8.2
    gd = Kd / (τd * s + 1)
    gp = K / (τ * s + 1)
    # P-control regulator first-order step response
    gc = Kc
    y_ovr_d = gd / (1 + gc * gp)
    @test isapprox(y_ovr_d, TransferFunction(Kd/(1+Kc*K) * [τ, 1], [τd*τ/(1+Kc*K), (τ+τd+Kc*K*τd)/(1+Kc*K), 1]))
    # PI-control servo first-order step response
    gc = Kc * (τI * s + 1) / (τI * s)
    y_ovr_y_sp = gc * gp / (1 + gc * gp)
    y_ovr_y_sp = pole_zero_cancellation(y_ovr_y_sp)
    @test isapprox(y_ovr_y_sp, TransferFunction([τI, 1], [τI*τ/(Kc*K), τI*(1+Kc*K)/(Kc*K), 1]), atol=1e-5)

    ###
    #  time delay
    ###
    g = exp(-4.0*s)
    @test g.time_delay == 4.0
    g = 1 / (s + 1) * exp(-6.2 * s)
    @test isapprox(g, TransferFunction([1], [1, 1], 6.2))
    @test isapprox(exp(-3*s), exp(-3.0*s))
    @test isapprox(exp(-s), TransferFunction([1], [1], 1.0))

    ###
    # zeros, poles, gain
    ###
    g = TransferFunction([3], [4, 3, 4, 1])
    z, p, k = zeros_poles_gain(g)
    @test k == 3
    g = TransferFunction([3], [4, 3, 4, 2])
    z, p, k = zeros_poles_gain(g)
    @test k == 3/2
    # G(s) = s + 2 / (s^2 - 4)
    tf = TransferFunction([1, 2], [1, 0, -4])
    z, p, k = zeros_poles_gain(tf)
    @test isapprox(k, -1/2, atol=1e-6)
    @test z == [-2.0]
    @test isapprox(p, [-2.0, 2.0])

    ###
    #  zeros, poles, k
    ###
    z = [-3.2, 4.0]
    p = [-3.0, -3.0, 0.0, 2.0]
    k = 23.2
    g = zeros_poles_k(z, p, k)
    _z, _p, _k = zeros_poles_k(g)
    sort!(_z)
    sort!(_p)
    @test isapprox(z, _z)
    @test isapprox(p, _p)
    @test isapprox(k, _k)
    g2 = zeros_poles_k(z, p, k, time_delay=3.0)
    @test ! isapprox(g, g2)
    @test isapprox(g2, g * exp(-3*s))
    g = 45.0 * s * s / (s * (s-23.0) * (s-8.0))
    _z, _p, _k = zeros_poles_k(g)
    sort!(_z)
    sort!(_p)
    @test isapprox(_k, 45.0)
    @test isapprox(_z, [0.0])
    @test isapprox(_p, [8.0, 23.0])
    p = [-2+im, -2-im] # try with complex numbers.
    z = [-1/5]
    @test isapprox(zeros_poles_k(z, p, 5.0), (5*s+1) / (s^2+4*s+5))
    @test ! isapprox(zeros_poles_k(z, p, 5.0, time_delay=1.0), (5*s+1) / (s^2+4*s+5))
    @test isapprox(zeros_poles_k(z, p, 5.0, time_delay=1.0), (5*s+1) / (s^2+4*s+5) * exp(-s))
    g = (6s^2+18*s+12)/(2*s^3+10*s^2+16*s+12)
    z, p, k = zeros_poles_k(g)
    sort!(z)
    @test isapprox(z, [-2.0, -1.0])
    @test isapprox(p, [-3.0, -1.0-im, -1+im])

    ###
    #  evaluate
    ###
    g = TransferFunction([1], [3, 1]) # 1 / (3s+1)
    z, p, k = zeros_poles_gain(g)
    @test evaluate(g, 0.0) == k
    @test evaluate(g, 1.0) == 1/4
    @test evaluate(g, 0.0) == 1.0
    # https://www.mathworks.com/help/control/ref/evalfr.html
    g = (s - 1) / (s ^ 2 + s + 1)
    z, p, k = zeros_poles_gain(g)
    @test isapprox(k, -1, atol=1e-6)
    @test evaluate(g, 0.0) == k
    @test isapprox(evaluate(g, 1+im), 0.2308 + 0.1538im, atol=0.0001)
    g = 1.0 / (s^2 + 2*s + 1)
    @test isapprox(evaluate(g, 0.1*im), 0.9705 - 0.1961im, atol=0.0001)

    ###
    #  properness
    ###
    tf = 1 / (2 * s + 1)
    @test proper(tf)
    @test strictly_proper(tf)
    tf = 3 * s / (2 * s + 1)
    @test proper(tf)
    @test ! strictly_proper(tf)
    tf = (3 * s * s) / (2 * s + 1)
    @test ! proper(tf)
    @test ! strictly_proper(tf)

    ###
    #  zeros, poles cancellation
    ###
    @test isapprox(pole_zero_cancellation((s-1)^2 / (s-1)), s-1)
    @test isapprox(pole_zero_cancellation((s-1) / (s-1)^2), 1/(s-1))
    @test isapprox(pole_zero_cancellation((2*s-2.0) / (s-1)^2), 2/(s-1))
    @test isapprox(pole_zero_cancellation(s * (9*s+3) / (s * (3*s-3))), (9*s+3)/(3*s-3), atol=1e-5)
    @test isapprox(pole_zero_cancellation(s^5/s), s^4)
    @test isapprox(pole_zero_cancellation(s^5*(s-1)/s/(s-1)), s^4)
    @test isapprox(pole_zero_cancellation(s^5*(s-1)/s/(s-1)*(s+2)), (s+2)*s^4)
    g = s * (s+1) / ((s+3) * s * (s+1) ^ 2)
    @test isapprox(pole_zero_cancellation(g, digits=4), 1 / ((s+3) * (s+1)))
    g = (s+1) / ((s^2-2*s+2) * (s+1) * s) # some imaginary poles
    @test isapprox(pole_zero_cancellation(g), 1/(s^2-2*s+2)/s)

    ###
    #  characteristic eqn.
    ###
    g_ol = 4 / (s + 3) / (s + 2) / (s + 1)
    @test isapprox(characteristic_polynomial(g_ol), Polynomial([10.0, 11, 6, 1], :s))

    ###
    #   order
    ###
    g = (s + 1) / (s + 1) ^ 2 # make sure it cancels the pole and zero.
    @test system_order(g) == (0, 1)
    @test system_order(first_order_system(1.0, 2.0)) == (0, 1)
    @test system_order(second_order_system(1.0, 2.0, 1.9)) == (0, 2)

    ###
    #  time constant, damping coefficient
    ###
    g = 4 / (6 * s + 2)
    @test isapprox(time_constant(g), 3.0)
    @test_throws ErrorException damping_coefficient(g)
    g = 1.0 / (8 * s^2 + 0.8 * s + 2)
    @test isapprox(time_constant(g), 2.0)
    @test_throws ErrorException time_constant(s / (s+4)) # only makes sense with first, second order...
    g = 1.0 / (8 * s^2 + 0.8 * s + 2)
    @test isapprox(damping_coefficient(g), 0.1) # 0.1
    τ = 4.3
    ξ = 1.2
    g = second_order_system(20.0, τ, ξ)
    @test isapprox(time_constant(g), τ)
    @test isapprox(damping_coefficient(g), ξ)
end

@testset "testing Simulation" begin
    # from http://web.mit.edu/2.14/www/Handouts/StateSpace.pdf
    tf = (13*s+26) / (1*s*s*s+7*s*s+19*s+13)
    A, B, C = Controlz.tf_to_ss(tf)
    @test isapprox(A, [0.0 1.0 0.0; 0  0 1; -13 -19 -7])
    @test isapprox(B, [0.0, 0.0, 1.0])
    @test isapprox(C, [26.0 13.0 0.0])

    # from https://www.engr.mun.ca/~millan/Eng6825/canonicals.pdf
    tf = TransferFunction([1, 3], [1, 3, 2])
    A, B, C = Controlz.tf_to_ss(tf)
    @test isapprox(A, [0.0 1.0; -2.0 -3.0])
    @test isapprox(B, [0.0, 1.0])
    @test isapprox(C, [3.0 1.0])

    # from https://people.kth.se/~demirel/State_Space_Representation_of_Transfer_Function_Systems.pdf
    tf = TransferFunction([1], [1, 6, 11, 6])
    A, B, C = Controlz.tf_to_ss(tf)
    @test isapprox(A, [0.0 1.0 0.0; 0.0 0.0 1.0; -6.0 -11.0 -6.0])
    @test isapprox(B, [0.0, 0.0, 1.0])
    @test isapprox(C, [1.0 0.0 0.0])

    tf = TransferFunction([1, 3, 3], [1, 2, 1])
    A, B, C = Controlz.tf_to_ss(tf)
    @test isapprox(A, [0.0 1.0; -1.0 -2.0])
    @test isapprox(B, [0.0, 1.0])
    @test isapprox(C, [2.0 1.0])

    # some inputs
    # L[cos(at)] = s/(s^2+a^2)
    a = 2.3
    U = s / (s^2+a^2)
    u_data = simulate(U, 12.0)
    _cosat(t::Float64) = t < 0.0 ? 0.0 : cos(a*t)
    @test isapprox(u_data[:, :output], _cosat.(u_data[:, :t]), rtol=0.001)
    # L[shifted step] = e^{-a*s}/s
    U = exp(-a*s) / s
    u_data = simulate(U, 12.0)
    _shifted_step(t::Float64) = t < a ? 0.0 : 1.0
    @test isapprox(u_data[:, :output], _shifted_step.(u_data[:, :t]), rtol=0.001)
    # L[t sin(at)] = 2as/(s^2+a^2)^2
    a = 2.1
    U = 2 * a * s / (s^2 + a^2)^2
    u_data = simulate(U, 12.0)
    u_truth = [t_i > 0.0 ? t_i * sin(a * t_i) : 0.0 for t_i in u_data[:, :t]]
    isapprox(u_data[:, :output], u_truth, rtol=0.01)
    # L{exp[-3(t-2)]S(t-2)} = e^{-2s} / (s+3)
    U = exp(-2*s) / (s + 3)
    u_data = simulate(U, 12.0)
    u_truth = [t_i > 2.0 ? exp(-3*(t_i - 2)) : 0.0 for t_i in u_data[:, :t]]
    isapprox(u_data[:, :output], u_truth, rtol=0.001)

    # first order step response
    K = 4.3
    τ = 2.8
    g = K / (τ * s + 1)
    u = 1 / s
    data = simulate(g * u, 12.0)
    y_truth = K * (1.0 .- exp.(-data[:, :t] ./ τ))
    y_truth[data[:, :t] .< 0.0] .= 0.0
    @test isapprox(y_truth, data[:, :output], rtol=0.0001)

    # first order impulse response
    data = simulate(g, 12.0)
    y_truth = K / τ * exp.(-data[:, :t]/τ)
    y_truth[data[:, :t] .< 0.0] .= 0.0
    @test isapprox(y_truth, data[:, :output], rtol=0.0001)

    # first order ramp input
    a = 2.0 # slope of ramp
    Y = g * a / s^2
    data = simulate(Y, 10.0)
    y_truth = K * a * data[:, :t] .+ K * a * τ * (exp.(- data[:, :t] ./ τ) .- 1.0)
    y_truth[data[:, :t] .< 0.0] .= 0.0
    @test isapprox(y_truth, data[:, :output], rtol=0.0001)

    # first order sinusoidal input
    ω = 2.0
    A = 4.5
    Y = g * A * ω / (s^2 + ω^2)
    data = simulate(Y, 10.0)
    ϕ = atan(-τ*ω)
    y_truth = τ * ω / (1 + (τ*ω)^2) * exp.(-data[:, :t] ./ τ) .+ 1 / sqrt(1+(τ*ω)^2) * sin.(ω*data[:, :t] .+ ϕ)
    y_truth *= K * A
    y_truth[data[:, :t] .< 0.0] .= 0.0
    @test isapprox(y_truth, data[:, :output], rtol=0.001)

    # FOPTD
    M = 3.3
    θ = 2.3
    τ = 0.7
    K = 9.0
    g = K / (τ * s + 1) * exp(-θ * s)
    u = M / s
    data = simulate(g * u, 10.0)
    y_truth = K * M * (1.0 .- exp.(-(data[:, :t] .- θ) ./ τ))
    y_truth[data[:, :t] .< θ] .= 0.0
    @test isapprox(y_truth, data[:, :output], atol=0.001)

    ##
    # second order, overdamped
    K = 4.3
    τ = 2.8
    ξ = 1.2
    g = K / (τ ^ 2 * s ^ 2 + 2 * τ * ξ * s + 1)
    # response to step
    M = 3.3
    Y = g * M / s
    data = simulate(Y, 20.0)
    y_truth = 1.0 .- exp.(- ξ / τ * data[:, :t]) .* (
        ξ / sqrt(ξ^2 - 1) * sinh.(sqrt.(ξ^2 - 1) / τ * data[:, :t]) .+ cosh.(sqrt.(ξ^2 - 1) / τ * data[:, :t]))
    y_truth *= K * M
    y_truth[data[:, :t] .< 0.0] .= 0.0
    @test isapprox(y_truth, data[:, :output], rtol=0.001)
    # response to impulse
    data = simulate(g, 20.0)
    y_truth = K / τ / sqrt(ξ^2-1) * exp.(-ξ/τ*data[:, :t]) .* sinh.(sqrt(ξ^2-1)/τ*data[:, :t])
    y_truth[data[:, :t] .< 0.0] .= 0.0
    @test isapprox(y_truth, data[:, :output], rtol=0.001)

    ##
    # second order, underdamped
    K = 4.3
    τ = 2.8
    ξ = 0.2
    g = K / (τ ^ 2 * s ^ 2 + 2 * τ * ξ * s + 1)
    # impulse response
    data = simulate(g, 20.0)
    y_truth = K / τ / sqrt(1-ξ^2) * exp.(-ξ/τ*data[:, :t]) .* sin.(sqrt(1-ξ^2)/τ*data[:, :t])
    y_truth[data[:, :t] .< 0.0] .= 0.0
    @test isapprox(y_truth, data[:, :output], rtol=0.001)

    ##
    # second order with zeros
    K = 32.0
    τₐ = 2.3
    τ₁ = 1.1
    τ₂ = 2.5
    g = K * (τₐ * s + 1) / (τ₁ * s + 1) / (τ₂ * s + 1)
    # step response
    data = simulate(g / s, 10.0) # direct inversion
    y_truth = 1 .- (τ₁ - τₐ) / (τ₁ - τ₂) * exp.(-data[:, :t]/τ₁) .- (τ₂ - τₐ) / (τ₂ - τ₁) * exp.(-data[:, :t]/τ₂)
    y_truth *= K
    y_truth[data[:, :t] .< 0.0] .= 0.0
    @test isapprox(y_truth, data[:, :output], rtol=0.001)

    # inerpolate
    t = [0.0, 1.0, 2.0]
    y = [0.0, 5.0, 10.0]
    data = DataFrame(t=t, output=y)
    @test isapprox(interpolate(data, 1.0), 5.0)
    @test isapprox(interpolate(data, 1.2), 1.2*5)
    τ = 3.45
    g = 1 / (τ * s + 1)
    data = simulate(g / s, 10.0)
    @test isapprox(interpolate(data, τ), 0.632, atol=0.002)
end

@testset "testing Controls" begin
    @test isapprox(TransferFunction(PIController(1.3, 2.0)), 1.3 * (1+1/(2*s)))
    @test isapprox(TransferFunction(PIDController(1.3, 2.0, 0.0)), 1.3 * (1+1/(2*s)))
    @test isapprox(TransferFunction(PIDController(1.3, 2.0, 0.1)), 1.3 * (1+1/(2*s) + 0.1*s))
end

@testset "testing margins" begin
    g_ol = 2 * exp(-s) / (5 * s + 1)
    m = gain_phase_margins(g_ol)
    @test isapprox(m.ω_c, 1.6887, atol=0.0001)
    @test isapprox(m.ω_g, 0.3464, atol=0.0001)
    @test isapprox(m.gain_margin, 4.2512, atol=0.0001)
    @test isapprox(m.phase_margin, 100.1535 / 180 * π, atol=0.0001)

    g_ol = TransferFunction([0.25], [1.0, 2, 1, 1])
    m = gain_phase_margins(g_ol)
    @test isapprox(m.ω_c, 1.0, atol=0.0001)
    @test isnan(m.ω_g)
    @test isapprox(m.gain_margin, 4.0, atol=0.0001)
    @test isnan(m.phase_margin)

    g_ol = 20*(1+1/(3*s)) * 4 /(s+4)/(s+6)*exp(-0.001*s)/(s+2)
    m = gain_phase_margins(g_ol)
    @test isapprox(m.ω_c, 6.32, atol=0.01)
    @test isapprox(m.gain_margin, 5.397, atol=0.01)

    # check asymptote behavior
    g_ol = 1 / (s + 1) # ω_c d.n.e.
    m = gain_phase_margins(g_ol)
    @test isnan(m.ω_c)
    @test isnan(m.gain_margin)
end

@testset "closed loop stuff" begin
    g_d =  3 / (s + 1) * exp(-5 * s)
    g_ol = 18 / (100 * s + 1) * exp(-s)
    cl = ClosedLoopTransferFunction(g_d, g_ol)

    # test conversion into standard form
    cls = Controlz.CLTFStandard(cl)
    @test cls.p_a ≈ Polynomial([3, 300], :s)
    @test cls.p_b ≈ Polynomial([1, 101, 100], :s)
    @test cls.p_c ≈ Polynomial([18, 18], :s)
    @test cls.ϕ ≈ 1.0
	@test cls.θ ≈ 5.0

    @test Controlz.order(cl) == 2
    @test Controlz.strictly_proper(cl)
    @test ! strictly_proper(ClosedLoopTransferFunction(g_d, g_ol * (4 * s + 1)))

    # algebra
    cl2 = cl * 4 / s
    @test cl2.g_ol == cl.g_ol
    @test cl2.top == cl.top * 4 / s

    cl3 = ClosedLoopTransferFunction(1 / (s + 1) * exp(-3 * s), 5 / (s + 2) * exp(-s))
    @test (-cl3).g_ol == 5 / (s + 2) * exp(-s)
    @test (-cl3).top == -1 / (s + 1) * exp(-3 * s)

    # run without a time delay to compare to other simulate function.
    gp = 3 / (s + 2)
    gc = TransferFunction(PIController(1.0, 3.0))
    data_old = simulate(gp * gc / (1 + gp * gc), 10.0)
    data_new = simulate(ClosedLoopTransferFunction(gp * gc, gp * gc), 10.0)
    @test isapprox(data_old[:, :output], data_new[:, :output], atol=0.0001)
end

@testset "example notebook (also to generate images for docs)" begin
    include(joinpath("..", "examples", "examples.jl"))
    @test true
end

max_stars_repo_name: SimonEnsemble/Controlz

max_issues_repo_name: SimonEnsemble/Controlz

max_forks_repo_name: SimonEnsemble/Controlz


text:
@benchgroup "centrality" begin
  @benchgroup "graphs" begin
    for (name, g) in GRAPHS
      @bench "$(name): degree" LightGraphs.degree_centrality($g)
      @bench "$(name): closeness" LightGraphs.closeness_centrality($g)
      if nv(g) < 1000
        @bench "$(name): betweenness" LightGraphs.betweenness_centrality($g)
        @bench "$(name): katz" LightGraphs.katz_centrality($g)
      end
    end
  end #graphs
  @benchgroup "digraphs" begin
    for (name, g) in DIGRAPHS
      @bench "$(name): degree" LightGraphs.degree_centrality($g)
      @bench "$(name): closeness" LightGraphs.closeness_centrality($g)
      if nv(g) < 1000
        @bench "$(name): betweenness" LightGraphs.betweenness_centrality($g)
        @bench "$(name): katz" LightGraphs.katz_centrality($g)
      end
      if nv(g) < 500
        @bench "$(name): pagerank"  LightGraphs.pagerank($g)
      end
    end
  end # digraphs
end # centrality

ext: jl


text:
module Mult_Test

using Compat
using Compat.Test
using IntervalArithmetic
using EAGOSmoothMcCormickGrad

function XaboutY(x,y,tol)
    return abs(x-y) <= tol
end

################################################################################
################################################################################
##############      Testing for Nonsmooth Standard Mult           ##############
################################################################################
################################################################################
EAGOSmoothMcCormickGrad.set_diff_relax(0)

################################################################################
################### Test Nonsmooth Zero in Both Case (Failing)   ###############
################################################################################
a = seed_g(Float64,1,2)
b = seed_g(Float64,2,2)
xIBox = [Interval(-2.0,1.0);Interval(-1.0,2.0)]
mBox = mid.(xIBox)
X = SMCg{2,Float64}(0.0,0.0,a,a,xIBox[1],false,xIBox,mBox)
Y = SMCg{2,Float64}(1.0,1.0,b,b,xIBox[2],false,xIBox,mBox)
out = X*Y

println("out: $out")

@test out.cc == 2.0
@test out.cv == -1.0
@test out.cc_grad[1] == 2.0
@test out.cc_grad[2] == -2.0
@test out.cv_grad[1] == 2.0
@test out.cv_grad[2] == 1.0
@test out.Intv.lo == -4.0
@test out.Intv.hi == 2.0

################################################################################
###################### Test Nonsmooth X1.l>0   (Passing)  ######################
################################################################################
a = seed_g(Float64,1,2)
b = seed_g(Float64,2,2)
xIBox = [Interval(1.0,5.0);Interval(-1.0,2.0)]
mBox = mid.(xIBox)
X = SMCg{2,Float64}(3.0,3.0,a,a,xIBox[1],false,xIBox,mBox)
Y = SMCg{2,Float64}(1.0,1.0,b,b,xIBox[2],false,xIBox,mBox)
out = X*Y

@test out.cc == 5.0
@test out.cv == 1.0
@test out.cc_grad[1] == 2.0
@test out.cc_grad[2] == 1.0
@test out.cv_grad[1] == 2.0
@test out.cv_grad[2] == 5.0
@test out.Intv.lo == -5.0
@test out.Intv.hi == 10.0

################################################################################
############## Test Nonsmooth X1.h<0  &&  X2.l>0 (Passing)  ######################
################################################################################

a = seed_g(Float64,1,2)
b = seed_g(Float64,2,2)
xIBox = [Interval(-6.0,-2.0);Interval(1.0,3.0)]
mBox = mid.(xIBox)
X = SMCg{2,Float64}(-4.0,-4.0,a,a,xIBox[1],false,xIBox,mBox)
Y = SMCg{2,Float64}(2.0,2.0,b,b,xIBox[2],false,xIBox,mBox)
out = X*Y

@test out.cc == -6.0
@test out.cv == -10.0
@test out.cc_grad[1] == 1.0
@test out.cc_grad[2] == -2.0
@test out.cv_grad[1] == 3.0
@test out.cv_grad[2] == -2.0
@test out.Intv.lo == -18.0
@test out.Intv.hi == -2.0

################################################################################
############## Test Nonsmooth X1.h<0  &&  X2.h<0 (Passing)  ######################
################################################################################

a = seed_g(Float64,1,2)
b = seed_g(Float64,2,2)
xIBox = [Interval(-6.0,-2.0);Interval(-7.0,-3.0)]
mBox = mid.(xIBox)
X = SMCg{2,Float64}(-4.0,-4.0,a,a,xIBox[1],false,xIBox,mBox)
Y = SMCg{2,Float64}(-5.0,-5.0,b,b,xIBox[2],false,xIBox,mBox)
out = X*Y

@test out.cc == 24.0
@test out.cv == 16.0
@test out.cc_grad[1] == -7.0
@test out.cc_grad[2] == -2.0
@test out.cv_grad[1] == -3.0
@test out.cv_grad[2] == -2.0
@test out.Intv.lo == 6.0
@test out.Intv.hi == 42.0

################################################################################
############## Test Nonsmooth X1.h<0  &&  0 in X2 (Passing)  ###################
################################################################################

a = seed_g(Float64,1,2)
b = seed_g(Float64,2,2)
xIBox = [Interval(-6.0,-2.0);Interval(-7.0,4.0)]
mBox = mid.(xIBox)
X = SMCg{2,Float64}(-4.0,-4.0,a,a,xIBox[1],false,xIBox,mBox)
Y = SMCg{2,Float64}(-5.0,-5.0,b,b,xIBox[2],false,xIBox,mBox)
out = X*Y

@test out.cc == 24.0
@test out.cv == 16.0
@test out.cc_grad[1] == -7.0
@test out.cc_grad[2] == -2.0
@test out.cv_grad[1] == -7.0
@test out.cv_grad[2] == -6.0
@test out.Intv.lo == -24.0
@test out.Intv.hi == 42.0

################################################################################
############## Test Nonsmooth 0 in X1  &&  X2.l > 0 ()  #################
################################################################################
a = seed_g(Float64,1,2)
b = seed_g(Float64,2,2)
xIBox = [Interval(-3.0,4.0);Interval(1.0,4.0)]
mBox = mid.(xIBox)
X = SMCg{2,Float64}(-2.0,-2.0,a,a,xIBox[1],false,xIBox,mBox)
Y = SMCg{2,Float64}(3.0,3.0,b,b,xIBox[2],false,xIBox,mBox)
out = X*Y

@test out.cc == -5.0
@test out.cv == -8.0
@test out.cc_grad[1] == 4.0
@test out.cc_grad[2] == -3.0
@test out.cv_grad[1] == 1.0
@test out.cv_grad[2] == -3.0
@test out.Intv.lo == -12.0
@test out.Intv.hi == 16.0

################################################################################
############## Test Nonsmooth 0 in X1  &&  X2.h < 0 ()         #################
################################################################################
a = seed_g(Float64,1,2)
b = seed_g(Float64,2,2)
xIBox = [Interval(-3.0,4.0);Interval(-5.0,-3.0)]
mBox = mid.(xIBox)
X = SMCg{2,Float64}(-2.0,-2.0,a,a,xIBox[1],false,xIBox,mBox)
Y = SMCg{2,Float64}(-4.0,-4.0,b,b,xIBox[2],false,xIBox,mBox)
out = X*Y

@test out.cc == 9.0
@test out.cv == 7.0
@test out.cc_grad[1] == -3.0
@test out.cc_grad[2] == -3.0
@test out.cv_grad[1] == -5.0
@test out.cv_grad[2] == -3.0
@test out.Intv.lo == -20.0
@test out.Intv.hi == 15.0

################################################################################
################################################################################
##############      Testing for Smooth Standard Mult           ##############
################################################################################
################################################################################
EAGOSmoothMcCormickGrad.set_diff_relax(1)

seed1 = seed_g(Float64,1,2)
seed2 = seed_g(Float64,2,2)
x1 = SMCg{2,Float64}(0.0,0.0,seed1,seed1,Interval(-200.0,200.0),false, [Interval(-200.0,200.0),Interval(0.0,400.0)],[0.0,200.0])
y1 = SMCg{2,Float64}(200.0,200.0,seed2,seed2,Interval(0.0,400.0),false, [Interval(-200.0,200.0),Interval(0.0,400.0)],[0.0,200.0])
z1 = x1*y1
@test XaboutY(z1.cc,40000,1E-4)
@test XaboutY(z1.cv,-40000,1E-4)

x2 = SMCg{2,Float64}(170.0,170.0,seed1,seed1,Interval(100.0,240.0),false, [Interval(100.0,240.0),Interval(100.0,400.0)],[170.0,250.0])
y2 = SMCg{2,Float64}(250.0,250.0,seed2,seed2,Interval(100.0,400.0),false, [Interval(100.0,240.0),Interval(100.0,400.0)],[170.0,250.0])
z2 = x2*y2
@test XaboutY(z2.cc,53000,1E-4)
@test XaboutY(z2.cv,32000,1E-4)

x3 = SMCg{2,Float64}(-200.0,-200.0,seed1,seed1,Interval(-300.0,-100.0),false, [Interval(-300.0,-100.0),Interval(-400.0,-200.0)],[-200.0,-300.0])
y3 = SMCg{2,Float64}(-300.0,-300.0,seed2,seed2,Interval(-400.0,-200.0),false, [Interval(-300.0,-100.0),Interval(-400.0,-200.0)],[-200.0,-300.0])
z3 = x3*y3
@test XaboutY(z3.cc,70000,1E-4)
@test XaboutY(z3.cv,50000,1E-4)

x4 = SMCg{2,Float64}(150.0,150.0,seed1,seed1,Interval(100.0,200.0),false, [Interval(100.0,200.0),Interval(-500.0,-100.0)],[150.0,-300.0])
y4 = SMCg{2,Float64}(-250.0,-250.0,seed2,seed2,Interval(-500.0,-100.0),false, [Interval(100.0,200.0),Interval(-500.0,-100.0)],[150.0,-300.0])
z4 = x4*y4
@test XaboutY(z4.cc,-30000,1E-3)
@test XaboutY(z4.cv,-47460.9375,1E-3)

x5 = SMCg{2,Float64}(-150.0,-150.0,seed1,seed1,Interval(-200.0,-100.0),false, [Interval(-200.0,-100.0),Interval(200.0,400.0)],[-150.0,300.0])
y5 = SMCg{2,Float64}(300.0,300.0,seed2,seed2,Interval(200.0,400.0),false, [Interval(-200.0,-100.0),Interval(200.0,400.0)],[-150.0,300.0])
z5 = x5*y5
@test XaboutY(z5.cc,-40000,1E-4)
@test XaboutY(z5.cv,-50000,1E-4)

end


text:
# fvalue.jl

#=
Metida.fvalue(lmm, [0 0 1 0 0 0; 0 0 0 1 0 0; 0 0 0 0 1 0])
=#
"""
    fvalue(lmm::LMM, l::Matrix)

F value for contrast matrix `l`.

```math
F = \\frac{\\beta'L'(LCL')^{-1}L\\beta}{rank(LCL')}
```
"""
function fvalue(lmm::LMM, l::Matrix)
    lcl = l*lmm.result.c*l'
    f   = lmm.result.beta'*l'*pinv(lcl)*l*lmm.result.beta/rank(lcl)
end


text:
using AtomicLevels

# Only capable of providing active sets for configurations with no
# holes below in any ell channel.
function active_set(ref_set::Config, i)
    open_set = filter(orb -> orb[4] == "*", ref_set)
    if i == 0
        return open_set
    end

    # Find all occupied n quantum numbers for each ell
    occ = Dict{Int,Vector{Int}}()
    for o in open_set
        ell = o[2]
        occ[ell] = push!(get(occ, ell, Vector{Int}([])), o[1])
    end

    ell_occ = keys(occ)
    ell_max = maximum(ell_occ)

    trunc_neg = v -> v > 0 ? v : 0

    active = Dict{Char,Int}()
    c = map(0:(ell_max+i)) do ell
        max_occ = maximum(get(occ, ell, [0]))
        # Amount of extra to add, if this ell channel is not occupied at all
        extra = trunc_neg(round(Int, ell>ell_max)*(i-(ell-ell_max)))
        n = max(max_occ + i, ell + extra + 1)
        Orbital((n,ell,1,"*"))
    end
end


text:
#=
  (Simple) True skill problem.

  Example inspired by
  Johannes Borgstrom, Andrew D. Gordon, Michael Greenberg, James Margetson, and Jurgen Van Gael:
  "Measure Transformer Semantics for Bayesian Machine Learning"
  https://www.microsoft.com/en-us/research/publication/measure-transformer-semantics-for-bayesian-machine-learning-2011/?from=http%3A%2F%2Fresearch.microsoft.com%2Fpubs%2F135344%2Fmsr-tr-2011-18.pdf

  Note that we constraint the generated worlds so they satisfies the
  observation constraints (a > b > c).

  Cf: ~/cplint/trueskill.pl
      ~/blog/true_skill.blog
      ~/psi/true_skill.psi
      ~/webppl/true_skill.wpl

   Compare with tug_of_war.jl

=#

#=
Summary Statistics
      parameters       mean       std   naive_se      mcse        ess      rhat
          Symbol    Float64   Float64    Float64   Float64    Float64   Float64

  performance[1]   104.1935    3.7222     0.0372    0.3157    42.6341    1.0651
  performance[2]   100.4037    2.3667     0.0237    0.1699   105.8531    0.9999
  performance[3]    96.8175    2.8633     0.0286    0.2260    58.5262    1.0032
        skill[1]   101.2267    2.6448     0.0264    0.1878    82.3399    1.0025
        skill[2]   100.5331    2.1493     0.0215    0.1536   124.2788    1.0023
        skill[3]    98.6982    2.5923     0.0259    0.2023    75.9888    1.0281

Quantiles
      parameters      2.5%      25.0%      50.0%      75.0%      97.5%
          Symbol   Float64    Float64    Float64    Float64    Float64

  performance[1]   98.8339   101.0689   103.6526   106.5544   111.5064
  performance[2]   95.9093    98.7352   100.5463   101.6910   104.9909
  performance[3]   92.2724    94.7470    96.3814    99.3117   102.2367
        skill[1]   95.0099    99.4469   101.5057   102.8886   105.4215
        skill[2]   96.3152    99.1755   100.4978   101.9342   104.8146
        skill[3]   93.1567    96.7794    99.0601   100.5969   103.7641

skill[1] > skill[2]
0.606
skill[1] > skill[3]
0.7464
skill[2] > skill[3]
0.6406
performance[1] > performance[2]
0.9997
performance[1] > performance[3]
1.0
performance[2] > performance[3]
1.0
  1.331539 seconds (7.62 M allocations: 386.380 MiB, 8.34% gc time)

=#

using Turing, StatsPlots, DataFrames
include("jl_utils.jl")

@model function true_skill(num_people)
    # There are three people, a, b, and c
    # num_people = 5

    # Each person has an unknown (latent) Skill and a
    # known performance, where the skill is
    # reflected in the performance (with uncertainties).
    # skill = Vector{Float64}(undef, num_people)
    # performance = Vector{Float64}(undef, num_people)
    skill = TArray{Float64}(undef, num_people)
    performance = TArray{Float64}(undef, num_people)
    for p in 1:num_people
        skill[p] ~ Normal(100,sqrt(10))
        performance[p] ~ Normal(skill[p], sqrt(15))
    end

    # Nicer (and slower), but it's not supported by IS or PG!
    # skill .~ Normal(100,sqrt(10))
    # @. performance ~ Normal(skill, sqrt(15)) # Note The @. since we refer to skill[p]

    for p1 in 1:num_people-1, p2 in p1+1:num_people
        true ~ Dirac(performance[p1] > performance[p2])
    end

    # What are their performance and (underlying) skills?

end

num_people = 5
model = true_skill(num_people)
# It seems that MH() is the only sampler that can handle this
# chains = sample(model, MH(), 10000) #
num_chains = 4
chains = sample(model, MH(), MCMCThreads(), 10_000, num_chains)

# chains = sample(model, PG(20), MCMCThreads(), 1000, num_chains) # PG doesn't support vectorizing assume statement
# chains = sample(model, SMC(1000), MCMCThreads(), 10_000, num_chains)

# Note: IS don't generate chains the same way as MH, PG, and SMC!
# chains = sample(model, IS(), MCMCThreads(), 1000, num_chains) # IS doesn't support vectorizing assume statement

# NUTS, HMC, and HMCDA throws this error:
# chains = sample(model, NUTS(), 1000) # Error
# chains = sample(model, HMC(0.1, 5), 1000) # Same error as NUTS
# chains = sample(model, HMCDA(0.15, 0.65), 1000) # Same errors as NUTS

display(chains)

# display(plot(chains))
# display(histogram(chains))
# display(gelmandiag(chains)) # requires number of chains >= 2

df = DataFrame(chains);

# Show the latent skills
for p1 in 1:num_people-1, p2 in p1+1:num_people
    println("skill[$p1] > skill[$p2]: ", mean(df[!,"skill[$p1]"] .>df[!,"skill[$p2]"]))
end

# Show the observed performance
for p1 in 1:num_people-1, p2 in p1+1:num_people
    println("performance[$p1] > performance[$p2]:", mean(df[!,"performance[$p1]"] .>df[!,"performance[$p2]"]))
end

max_stars_repo_name: Wikunia/hakank

max_issues_repo_name: Wikunia/hakank

max_forks_repo_name: Wikunia/hakank


text:
# By starting at the top of the triangle below and moving to adjacent numbers
# on the row below, the maximum total from top to bottom is 23.
#
# 3
# 7 4
# 2 4 6
# 8 5 9 3
#
# That is, 3 + 7 + 4 + 9 = 23.
#
# Find the maximum total from top to bottom in triangle.txt (right click and
# 'Save Link/Target As...'), a 15K text file containing a triangle with
# one-hundred rows.
#
# NOTE: This is a much more difficult version of Problem 18. It is not
# possible to try every route to solve this problem, as there are 299
# altogether! If you could check one trillion (1012) routes every second it
# would take over twenty billion years to check them all. There is an efficient
# algorithm to solve it. ;o)

using ProjectEulerSolutions

# Convert a string representation of a pyramid of numbers from a file into an
# integer representation.
function pyramid_from_file(path::String)::Vector{Vector{Integer}}
    pyramid = Vector{Vector{Integer}}()
    open(path) do file
        for row in eachline(file)
            if length(row) > 0
                cols = map(x->convert(Integer, parse(Int, x)), split(row, " "))
                push!(pyramid, cols)
            end
        end
    end
    return pyramid
end

p067 = Problems.Problem(max_path_through_pyramid)

input_path = "data/p067_triangle.txt"

Problems.benchmark(p067, pyramid_from_file(input_path))

max_stars_repo_name: gnujosh/julia-euler-project

max_issues_repo_name: gnujosh/julia-euler-project

max_forks_repo_name: gnujosh/julia-euler-project


text:
export @t_str, @tex_str, @c_str

const TEX_STR_PAT  = r"##([_\p{L}](?:[\p{L}\d_]*))"

"""
    tex_str(s)

E/ Allows use of \$ and \\ while also allowing interpolation via `##name`.
This is useful for labels that contain tex-like expressions that may be
generated in a sequence.
"""
macro tex_str(s)
    m = match(TEX_STR_PAT, s)
    m === nothing && return s
    v = Symbol(m.captures[1])
    esc(:(replace($s, $TEX_STR_PAT => string(eval($v)))))
end

@eval const $(Symbol("@t_str")) = $(Symbol("@tex_str"))
@eval const $(Symbol("@c_str")) = $(Symbol("@colorant_str"))

const BRACKETS_PAT = r"[\(\),\.]"

"""
    col2str(c; str)

I/ Return a GLE-compatible string representation of a color.
In the case of Markerstyle color with the relevant routine it is useful
to get a representation devoid of brackets, dots or dashes.
"""
function col2str(
            col::Colorant;
            str=false
        )::String

    crgba = convert(RGBA, col)
    r, g, b, α = fl2str.(Float64.((crgba.r, crgba.g, crgba.b, crgba.alpha)))
    s = "rgba($r,$g,$b,$α)"
    str || return s
    return replace(s, BRACKETS_PAT => "_")
end

"""
    fl2str(f; d)

I/ Round a float `f` to `d` digits after the decimals and return the string.
"""
function fl2str(
            f::Float64;
            d::Int=4
        )::String

    return string(round(f, digits=d))
end
fl2str(r::Real; kw...) = fl2str(Float64(r); kw...)

"""
    vec2str(v; sep)

I/ GLE-compatible string representation of a vector.
"""
vec2str(v::AV{String}; sep=" ") = join(("\"$vi\""  for vi in v), sep)
vec2str(v::AV{<:Real}; sep=" ") = join(fl2str.(v), sep)
dlist(rge::UnitRange;  sep=" ") = join(("d$i" for i in rge), sep)


text:
module HarmonieNamelist

export schema2namelist

tofortran(val::Bool) =  val ? ".TRUE."  :  ".FALSE."
tofortran(val::Number) = val
tofortran(val::AbstractString) = "'$val'"
tofortran(val::Nothing) = ""

function schema2namelist(io::IO, d::AbstractDict)
    println(io, "# $(d["description"])")
    for link in d["links"]
        println(io,"# $(link["rel"])")
        println(io,"# $(link["href"])")
        println(io,"")
    end
    println(io,"&$(d["title"])   ")
    for (key,property) in d["properties"]
        default     = tofortran(get(property,"default",""))
        description = property["description"]
        keydefault  = rpad("$key=$default",25)
        enum        = join(tofortran.(get(property,"enum",nothing)),", ")
        println(io, "  $keydefault  # $description [$enum]")
    end
    println(io, "/")
    println(io,"")
end

end # module


text:
module TestCUDA

import LinearInterpolations as LI
import AxisArrayHelpers as AH
using CUDA
using Test

@testset "TextureEx" begin
    @testset "1d" begin
        axes_from = (CuVector(Float32[1,2,3,3.1, 5]),)
        axes_to = (1f0:4f0,)
        vals_to = CuVector(Float32[10,20,30,40])
        data_to = (axes=axes_to, values=vals_to)
        res_tex = AH.pullback(identity, axes_from, data_to, executor=AH.TextureEx(), extrapolate=LI.Replicate())
        @test only(res_tex.axes) == only(axes_from)
        @test res_tex.values isa CUDA.CuVector
        @test Vector(res_tex.values) ≈ [10,20,30,31,40] rtol=1e-3
    end
    @testset "2d" begin
        axes_from = (CuVector(Float32[1,2,3,3.1, 4, 5]),)
        axes_to = (-1f0:2:1f0, 1f0:4f0)
        vals_to = CuMatrix(
            Float32[
             10  20  30  40;
            100 200 300 400
        ])
        function f(pt)
            x, = pt
            if isodd(round(Int, x))
                (-1f0, x)
            else
                (1f0, x)
            end
        end

        data_to = (axes=axes_to, values=vals_to)
        res_tex = AH.pullback(f, axes_from, data_to, executor=AH.TextureEx(), extrapolate=LI.Replicate())
        @test Tuple(res_tex.axes) == axes_from
        @test res_tex.values isa CUDA.CuVector
        @test Vector(res_tex.values) ≈ [10,200,30,31,400, 40] rtol=1e-3
    end

    @testset "3d" begin
        function f(pt)
            x,y,z=pt
            z+sin(x), y+cos(z), x+1f-3*sqrt(abs(y*z))
        end
        for _ in 1:10
            axes_from = (-10f0:2:10f0, -10f0:1:12f0, -12f0:3:9f0)
            axes_to = (-8:1f0:8, -9:2f0:11f0, -20:1f0:10)
            vals_to = randn(Float32, map(length, axes_to))
            data_to = (axes=axes_to, values=vals_to)
            res_cpu = AH.pullback(f, axes_from, data_to, executor=AH.ThreadsEx(), extrapolate=LI.Replicate())
            res_tex = AH.pullback(f, axes_from, data_to, executor=AH.TextureEx(), extrapolate=LI.Replicate())
            @test res_cpu.axes === res_tex.axes
            @test res_tex.values isa CUDA.CuArray
            @test res_cpu.values ≈ Array(res_tex.values) rtol=5e-3
        end
    end
end

end#module


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

mutable struct ApiError <: SwaggerModel
    details::Any # spec type: Union{ Nothing, Vector{ApiErrorBase} } # spec name: details
    innererror::Any # spec type: Union{ Nothing, InnerError } # spec name: innererror
    code::Any # spec type: Union{ Nothing, String } # spec name: code
    target::Any # spec type: Union{ Nothing, String } # spec name: target
    message::Any # spec type: Union{ Nothing, String } # spec name: message

    function ApiError(;details=nothing, innererror=nothing, code=nothing, target=nothing, message=nothing)
        o = new()
        validate_property(ApiError, Symbol("details"), details)
        setfield!(o, Symbol("details"), details)
        validate_property(ApiError, Symbol("innererror"), innererror)
        setfield!(o, Symbol("innererror"), innererror)
        validate_property(ApiError, Symbol("code"), code)
        setfield!(o, Symbol("code"), code)
        validate_property(ApiError, Symbol("target"), target)
        setfield!(o, Symbol("target"), target)
        validate_property(ApiError, Symbol("message"), message)
        setfield!(o, Symbol("message"), message)
        o
    end
end # type ApiError

const _property_map_ApiError = Dict{Symbol,Symbol}(Symbol("details")=>Symbol("details"), Symbol("innererror")=>Symbol("innererror"), Symbol("code")=>Symbol("code"), Symbol("target")=>Symbol("target"), Symbol("message")=>Symbol("message"))
const _property_types_ApiError = Dict{Symbol,String}(Symbol("details")=>"Vector{ApiErrorBase}", Symbol("innererror")=>"InnerError", Symbol("code")=>"String", Symbol("target")=>"String", Symbol("message")=>"String")
Base.propertynames(::Type{ ApiError }) = collect(keys(_property_map_ApiError))
Swagger.property_type(::Type{ ApiError }, name::Symbol) = Union{Nothing,eval(Meta.parse(_property_types_ApiError[name]))}
Swagger.field_name(::Type{ ApiError }, property_name::Symbol) =  _property_map_ApiError[property_name]

function check_required(o::ApiError)
    true
end

function validate_property(::Type{ ApiError }, name::Symbol, val)
end


text:
# MIT License
# Copyright (c) 2019: Antonio Saragga Seabra

# 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.

"""
    nfcaixa(startdate::Date, maturity::Date, frequency)
        startdate: date scalar
        maturity: date scalar
        frequency: 0 (zero-coupn), 1 (annual), 2 , 3, 4, 6, 12 (monthly).

Computes the **number of cash-flows** of a bond between dates [startdate, maturity], given the annual frequency - freq.
"""
function nfcaixa(startdate::Date, enddate::Date, frequency::Int64)
    if freq == 0
        return 1
    elseif rem(12,frequency) != 0
        return println("\x1b[94m*** ERROR: freq = $freq. Note that freq must be: 0 (zero-coupon), 1 (annual), 2, 3, 4, 6, 12 (monthly).\x1b[37m")
    else
        i = 0
        date = startdate
        while date <= enddate
            i = i + 1
            date = date + Dates.Month(12/frequency)
        end
        return i
    end
end

"""
Bond Cashflows - **Bullet Bonds**

function cf(startdate::Date, maturity::Date, couponrate, frequency, [principal = 100])
    startdate: date scalar
    maturity: date scalar
    frequency: 0 (zero-coupn), 1 (annual), 2 , 3, 4, 6, 12 (monthly).

    kx2 array - Calculates the timing and value (interest + principal) of all cash flows. The entire principal value is paid all at once on the maturity date. Principal is by default set equal to 100. Zero-coupons have 'freq' equal to 0.
"""
function fcaixa(startdate::Date, enddate::Date, couponrate, frequency::Int64; principal = 100)
    if frequency == 0
        return [enddate principal]
    elseif rem(12,frequency) != 0
        return println("\x1b[94m*** ERROR: freq = $freq. Note that freq must be: 0 (zero-coupon), 1 (annual), 2, 3, 4, 6, 12 (monthly).\x1b[37m")
    else
        cfvalue = (1+couponrate/frequency) * principal
        cashflows = [enddate cfvalue]
        cfdate = enddate - Dates.Month(12/frequency)
        while startdate < cfdate
            cashflows = [cashflows; [cfdate couponrate/frequency*principal]]
            cfdate = cfdate - Dates.Month(12/frequency)
        end
        return cashflows
    end
end

"""
Bond Cashflows - **Bonds with Amortization of Principal**

function fcaixa(startdate::Date, maturity::Date, couponrate, frequency, amortiz, [principal = 100])
    startdate: date scalar
    maturity: date scalar
    frequency: 0 (zero-coupn), 1 (annual), 2 , 3, 4, 6, 12 (monthly).
    amortiz: scalar

    kx2 Array{Date,Float64} - Calculates the timing and value (interest+principal) of all cash flows. The entire principal value is paid all at once on the maturity date. There is a linear amortization schedulle.
"""
function fcaixa(startdate::Date, enddate::Date, couponrate, frequency::Int64, amortiz; principal = 100)
    if rem(principal,amortiz) != 0
        println("\x1b[94m*** ERROR: Principal is not an integer multiple of the amortization value!\x1b[37m")
    end
    i = 1
    cfvalue = (1+couponrate/frequency) * amortiz
    cashflows = [enddate cfvalue]
    cfdate = enddate - Dates.Month(12/frequency)
    while startdate < cfdate
        i = i + 1
        cashflows = [cashflows; [cfdate (1 + couponrate/frequency * i ) * amortiz]]
        cfdate = cfdate - Dates.Month(12/frequency)
    end
    return cashflows
end

"""
spreço(cf, r, m, dc)

Computes the clean price of a bond
    r:  scalar or K vector of interest rates
    cf: scalar or K vector of cash flows
    m: K vector of times for the cash flows
    dc: Discrete/Continuous composition
"""
function spreço
end
function spreço(cf,m,r, ::Type{Discrete})  #cf is a vector of all cash flows at times m
    dcf = cf./((r.+1).^m)                     #cf1/(1+r1)^m1 + cf2/(1+r2)^m2 + ...
    return sum(dcf)
end

function spreço(cf,m,r, ::Type{Continuous})  #cf is a vector of all cash flows at times m
    dcf = cf.* exp.(-r.*m)                      #cf1*exp(-r1*m1) + cf2*exp(-r2*m2) + ...
    return sum(dcf)
end

"""
spreço(g, k, h, n, f1, f2, C, v, y)

Computes the dirty price (clean price plus accrued interest) of the bond using the ISMA formula
Brown, P. J. (1998): Bond Markets, Structures and Yields Calculations. ISMA Publication.
"""
function spreço(g, k, h, n, f1, f2, C, v, y)
    v = 1/(1+y/h)
    return v^f1 * (k + g/h * (v * (1-v^(n-1)) / (1-v)) + (C + g/h * f_2)*v^(n+f2-1))
end

"""
spreço(cf, m, r, z, dc)

Computes the dirty price (clean price plus accrued interest) of the bond using information about the z-spread and benchmark interest rates
    r:  scalar or K vector of benchmark interest rates
    cf: scalar or K vector of cash flows
    m:  K vector of times for the cash flows
    z:  scalar, the z-spread
    dc: Discrete/Continuous composition
"""
function spreço(cf,m,r, z, ::Type{Discrete})  #cf is a vector of all cash flows at times m
    dcf = cf./((r.+1).^m .*(z+1).^m)                     #cf1/(1+r1)^m1 + cf2/(1+r2)^m2 + ...
    return sum(dcf)
end
function spreço(cf,m,r, z, ::Type{Continuous})  #cf is a vector of all cash flows at times m
    dcf = cf.* exp.(-(r.+z).*m)                      #cf1*exp(-r1*m1) + cf2*exp(-r2*m2) + ...
    return sum(dcf)
end


text:
#Método de Euler
function Euler(df, f0, ti, tf, n)
    f = Array{Float64}(n)
    h = (tf - ti) / n
    f[1] = f0 + h*df(f0, ti)
    tk = ti + h
    for k = 1:(n-1)
        f[k + 1] = f[k] + h*df(f[k], tk)
        tk = tk + h
    end
    return vcat(f0, f)
end

#Método de Euler Modificado ou RUnge Kutta 2ª Orde
function Runge_Kutta_2(df, f0, ti, tf, n)
    f = Array{Float64}(n)
    h = (tf - ti) / n
    k1 = df(f0, ti)
    k2 = df(f0 + k1*h/2, ti + h/2)
    f[1] = f0 + k2*h
    tk = ti + h
    for k = 1:(n-1)
        k1 = df(f[k], tk)
        k2 = df(f[k] + k1*h/2, tk + h/2)
        f[k+1] = f[k] + k2*h
        tk = tk + h
    end
    return vcat(f0, f)
end

#Método de RUnge Kutta 3ª Ordem
function Runge_Kutta_3(df, f0, ti, tf, n)
    f = Array{Float64}(n)
    h = (tf - ti) / n
    k1 = h*df(f0, ti)
    k2 = h*df(f0 + k1/2, ti + h/2)
    k3 = h*df(f0 + k2*3/4, ti + 3*h/4)
    f[1] = f0 + (2/9)*k1 + (1/3)*k2 + (4/9)*k3
    tk = ti + h
    for k = 1:(n-1)
        k1 = h*df(f[k], tk)
        k2 = h*df(f[k] + k1/2, tk + h/2)
        k3 = h*df(f[k] + k2*3/4, tk + h/2)
        f[k + 1] = f[k] + (2/9)*k1 + (1/3)*k2 + (4/9)*k3
        tk = tk + h
    end
    return vcat(f0, f)
end

#Método de RUnge Kutta 4ª Ordem
function Runge_Kutta_4(df, f0, ti, tf, n)
    f = Array{Float64}(n)
    h = (tf - ti) / n
    k1 = h*df(f0, ti)
    k2 = h*df(f0 + k1/2, ti + h/2)
    k3 = h*df(f0 + k2/2, ti + h/2)
    k4 = h*df(f0 + k3, ti + h)
    f[1] = f0 + (k1 + 2*(k2 + k3) + k4)*(1/6)
    tk = ti + h
    for k = 1:(n-1)
        k1 = h*df(f[k], tk)
        k2 = h*df(f[k] + k1/2, tk + h/2)
        k3 = h*df(f[k] + k2/2, tk + h/2)
        k4 = h*df(f[k] + k3, tk + h)
        f[k + 1] = f[k] + (k1 + 2*(k2 + k3) + k4)*(1/6)
        tk = tk + h
    end
    return vcat(f0, f)
end

using Plots

#Equação Diferencial Ordinária
df(y, t) = -t/y
f0 = 20
ti = 0
tf = 16
n = 4

#Comparação dos Métodos Com a Resposta Analítica
y = collect(ti:(tf-ti)/n:tf)
y1(y) = 0.1*exp.(2*y)

y2 = Euler(df, f0, ti, tf, n)
y3 = Runge_Kutta_2(df, f0, ti, tf, n)
y4 = Runge_Kutta_3(df, f0, ti, tf, n)
y5 = Runge_Kutta_4(df, f0, ti, tf, n)
plot(y, y1,  label="Analítico")
plot!(y, y2, label="Euler")
plot!(y, y3, label="Runge Kutta 2")
plot!(y, y4, line=(:dot,3), label="Runge Kutta 3")
plot!(y, y5, line=(:dot,4), label="Runge Kutta 4")

#Análise do Erro Relativo
erro1 = abs(y1(y) - y2)
erro2 = abs(y1(y) - y3)
erro3 = abs(y1(y) - y4)
erro4 = abs(y1(y) - y5)
plot(y, erro1, label="Erro Euler")
plot!(y, erro2, label="Erro Runge Kutta 2")
plot!(y, erro3, line=(:dot,3), label="Erro Runge Kutta 3")
plot!(y, erro4, line=(:dot,4), label="Ero Runge Kutta 4")

max_stars_repo_path: 001 - Equacao_Diferencia_Com_PVI.jl

max_stars_repo_name: LudwigFonseca/Julia-Language-Examples

max_issues_repo_path: 001 - Equacao_Diferencia_Com_PVI.jl

max_issues_repo_name: LudwigFonseca/Julia-Language-Examples

max_forks_repo_path: 001 - Equacao_Diferencia_Com_PVI.jl

max_forks_repo_name: LudwigFonseca/Julia-Language-Examples


text:
# This file includes macro tools to define new components types

function _append_common_fields!(ex, newbody, newparamtypes)
    # Append body
    body = ex.args[3]
    append!(body.args, newbody.args)

    # Append struct type parameters
    name = ex.args[2]
    if name isa Expr && name.head === :(<:)
        name = name.args[1]
    end

    if name isa Expr && name.head === :curly
        append!(name.args, newparamtypes)
    elseif name isa Symbol
        ex.args[2] = Expr(:curly, name, newparamtypes...)  # parametrize ex
    end
end

function deftype(ex)
    # Check ex head
    ex isa Expr && ex.head == :struct || error("Invalid source defition")

    # Get struct name
    name = ex.args[2]
    if name isa Expr && name.head === :(<:)
        name = name.args[1]
    end

    # Process struct body
    body = ex.args[3]
    kwargs = Expr(:parameters)
    callargs = Symbol[]

    _def!(body, kwargs, callargs)

    # struct has no fields
    isempty(kwargs.args) && return quote
        Base.@__doc__($(esc(ex)))
    end

    if name isa Symbol
        return quote
            Base.@__doc__($(esc(ex)))
            $(esc(name))($kwargs) = $(esc(name))($(callargs...))
        end
    elseif name isa Expr && name.head === :curly
        _name = name.args[1]
        _param_types = name.args[2:end]
        __param_types = [_type_ isa Symbol ? _type_  : _type_.args[1] for _type_ in _param_types]
        return quote
            Base.@__doc__($(esc(ex)))
            $(esc(_name))($kwargs) = $(esc(_name))($(callargs...))
            $(esc(_name)){$(esc.(__param_types)...)}($kwargs) where {$(esc.(_param_types)...)} =
                $(esc(_name)){$(esc.(__param_types)...)}($(callargs...))
        end
    end
end

function _def!(body, kwargs, callargs)
    for i in 1 : length(body.args)
        bodyex = body.args[i]
        if bodyex isa Symbol # var
            push!(kwargs.args, bodyex)
            push!(callargs, bodyex)
        elseif bodyex isa Expr
            if bodyex.head === :(=)
                rhs = bodyex.args[2]
                lhs = bodyex.args[1]
                if lhs isa Expr && lhs.head === :(::) # var::T = 1
                    var = lhs.args[1]
                elseif lhs isa Symbol # var = 1
                    var = lhs
                elseif lhs isa Expr && lhs.head == :call # inner constructors
                    continue
                end
                push!(kwargs.args, Expr(:kw, var, esc(rhs)))
                push!(callargs, var)
                body.args[i] = lhs
            elseif bodyex.head === :(::)  # var::T
                var = bodyex.args[1]
                push!(kwargs.args, var)
                push!(callargs, var)
            end
        end
    end
end


text:
# TODO: find/replace all [PkgName] with CamelCase, all [pkgname] with lowercase

# [WEBSITE]

function _initialize_backend(::[PkgName]AbstractBackend; kw...)
  @eval begin
    import [PkgName]
    export [PkgName]
    # TODO: other initialization that needs to be eval-ed
  end
  # TODO: other initialization
end

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

function _create_plot(pkg::[PkgName]AbstractBackend, d::KW)
  # TODO: create the window/canvas/context that is the plot within the backend (call it `o`)
  # TODO: initialize the plot... title, xlabel, bgcolor, etc
  Plot(nothing, pkg, 0, d, KW[])
end

function _add_series(::[PkgName]AbstractBackend, plt::Plot, d::KW)
  # TODO: add one series to the underlying package
  push!(plt.seriesargs, d)
  plt
end

function _add_annotations{X,Y,V}(plt::Plot{[PkgName]AbstractBackend}, anns::AVec{@compat(Tuple{X,Y,V})})
  for ann in anns
    # TODO: add the annotation to the plot
  end
end

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

function _before_update_plot(plt::Plot{[PkgName]AbstractBackend})
end

# TODO: override this to update plot items (title, xlabel, etc) after creation
function _update_plot(plt::Plot{[PkgName]AbstractBackend}, d::KW)
end

function _update_plot_pos_size(plt::AbstractPlot{[PkgName]AbstractBackend}, d::KW)
end

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

# accessors for x/y data

# function getxy(plt::Plot{[PkgName]AbstractBackend}, i::Int)
#   # TODO: return a tuple of (x, y) vectors
# end
#
# function setxy!{X,Y}(plt::Plot{[PkgName]AbstractBackend}, xy::Tuple{X,Y}, i::Integer)
#   # TODO: set the plot data from the (x,y) tuple
#   plt
# end

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

function _create_subplot(subplt::Subplot{[PkgName]AbstractBackend}, isbefore::Bool)
  # TODO: build the underlying Subplot object.  this is where you might layout the panes within a GUI window, for example
end

function _expand_limits(lims, plt::Plot{[PkgName]AbstractBackend}, isx::Bool)
  # TODO: call expand limits for each plot data
end

function _remove_axis(plt::Plot{[PkgName]AbstractBackend}, isx::Bool)
  # TODO: if plot is inner subplot, might need to remove ticks or axis labels
end

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

function Base.writemime(io::IO, ::MIME"image/png", plt::AbstractPlot{[PkgName]AbstractBackend})
  # TODO: write a png to io
end

function Base.display(::PlotsDisplay, plt::Plot{[PkgName]AbstractBackend})
  # TODO: display/show the plot
end

function Base.display(::PlotsDisplay, plt::Subplot{[PkgName]AbstractBackend})
  # TODO: display/show the subplot
end


text:
module AprilTags

using AprilTags_jll
using Requires
using DocStringExtensions
using LinearAlgebra, Statistics
using Colors, ImageDraw, FixedPointNumbers
import Base.convert

export
#helpers
AprilTag,
AprilTagDetector,
freeDetector!,
getTagDetections,
homography_to_pose,
homographytopose,
threadcalldetect,
getAprilTagImage,
detectAndPose,
tagOrthogonalIteration,
# wrappers
apriltag_detector_create,
tag36h11_create,
tag36h11_destroy,

apriltag_detector_add_family,
apriltag_detector_detect,
apriltag_detections_destroy,
apriltag_detector_destroy,
threadcall_apriltag_detector_detect,

#drawing and plotting
drawTagBox!,
drawTagAxes!

include("wrapper.jl")
include("helpers.jl")
include("tagdraw.jl")
include("additionalutils.jl")
include("calibrationutils.jl")

function __init__()
    # conditional requirement
    @require FreeTypeAbstraction="663a7486-cb36-511b-a19d-713bb74d65c9" include("tagtext.jl")
end

end # module