text:
@testset "Field types" begin
    @testset "Linear fields" begin
        @field(F) do
            I₀ = 1.0
            T = 2.0
            σ = 3.0
            Tmax = 3.0
        end

        @test parent(F) === F

        @test field_envelope(F, 0.0) ≈ 1 rtol=1e-7
        @test fluence(F) ≈ 5.452488660809239*3.0*√(2π)/π

        withenv("UNITFUL_FANCY_EXPONENTS" => true) do
            @test string(F) == """
                               Linearly polarized field with
                                 - I₀ = 1.0000e+00 au = 3.5094452e16 W cm⁻² =>
                                   - E₀ = 1.0000e+00 au = 514.2207 GV m⁻¹
                                   - A₀ = 0.3183 au
                                 – a Fixed carrier @ λ = 14.5033 nm (T = 48.3777 as, ω = 3.1416 Ha = 85.4871 eV)
                                 – and a Gaussian envelope of duration 170.8811 as (intensity FWHM; ±2.00σ)
                                 – Uₚ = 0.0253 Ha = 689.2724 meV => α = 0.1013 Bohr = 5.3617 pm"""
        end
        withenv("UNITFUL_FANCY_EXPONENTS" => false) do
            @test string(F) == """
                               Linearly polarized field with
                                 - I₀ = 1.0000e+00 au = 3.5094452e16 W cm^-2 =>
                                   - E₀ = 1.0000e+00 au = 514.2207 GV m^-1
                                   - A₀ = 0.3183 au
                                 – a Fixed carrier @ λ = 14.5033 nm (T = 48.3777 as, ω = 3.1416 Ha = 85.4871 eV)
                                 – and a Gaussian envelope of duration 170.8811 as (intensity FWHM; ±2.00σ)
                                 – Uₚ = 0.0253 Ha = 689.2724 meV => α = 0.1013 Bohr = 5.3617 pm"""
        end

        @test ElectricFields.rotation_matrix(F) == Matrix(I, 3, 3)

        t = [0.4, 0.5]
        for fun in (vector_potential, field_amplitude, intensity)
            @test fun(F, t) == [fun(F, t[1]), fun(F, t[2])]
        end
    end

    @testset "Transverse fields" begin
        @field(F) do
            I₀ = 2.0
            T = 2.0
            σ = 3.0
            Tmax = 3.0
            rotation = π/4, [0,0,1]
        end

        @test parent(F) === F

        @test field_envelope(F, 0.0) ≈ √2 rtol=1e-7
        @test fluence(F) ≈ 5.452488660809239*2*3.0*√(2π)/π

        @test intensity(F) == 2.0
        @test amplitude(F) ≈ √2 rtol=1e-7

        withenv("UNITFUL_FANCY_EXPONENTS" => true) do
            @test string(F) == """
                               Transversely polarized field with
                                 - I₀ = 2.0000e+00 au = 7.0188904e16 W cm⁻² =>
                                   - E₀ = 1.4142e+00 au = 727.2178 GV m⁻¹
                                   - A₀ = 0.4502 au
                                 – a LinearTransverseCarrier: Fixed carrier @ λ = 14.5033 nm (T = 48.3777 as, ω = 3.1416 Ha = 85.4871 eV)
                                 – a Gaussian envelope of duration 170.8811 as (intensity FWHM; ±2.00σ)
                                 – and a rotation of 0.25π about [0.000, 0.000, 1.000]
                                 – Uₚ = 0.0507 Ha = 1.3785 eV => α = 0.1433 Bohr = 7.5826 pm"""
        end
        withenv("UNITFUL_FANCY_EXPONENTS" => false) do
            @test string(F) == """
                               Transversely polarized field with
                                 - I₀ = 2.0000e+00 au = 7.0188904e16 W cm^-2 =>
                                   - E₀ = 1.4142e+00 au = 727.2178 GV m^-1
                                   - A₀ = 0.4502 au
                                 – a LinearTransverseCarrier: Fixed carrier @ λ = 14.5033 nm (T = 48.3777 as, ω = 3.1416 Ha = 85.4871 eV)
                                 – a Gaussian envelope of duration 170.8811 as (intensity FWHM; ±2.00σ)
                                 – and a rotation of 0.25π about [0.000, 0.000, 1.000]
                                 – Uₚ = 0.0507 Ha = 1.3785 eV => α = 0.1433 Bohr = 7.5826 pm"""
        end

        @test ElectricFields.rotation_matrix(F) ≈ [1/√2 -1/√2 0
                                                   1/√2 1/√2  0
                                                   0    0     1]

        t = [0.4, 0.5]
        for fun in (vector_potential, field_amplitude, intensity)
            @test fun(F, t) == [fun(F, 0.4)'; fun(F, 0.5)']
        end

        @field(F2) do
            I₀ = 2.0
            T = 2.0
            σ = 3.0
            Tmax = 3.0
            ξ = 1.0
        end

        withenv("UNITFUL_FANCY_EXPONENTS" => true) do
            @test string(F2) == """
                                Transversely polarized field with
                                  - I₀ = 2.0000e+00 au = 7.0188904e16 W cm⁻² =>
                                    - E₀ = 1.4142e+00 au = 727.2178 GV m⁻¹
                                    - A₀ = 0.4502 au
                                  – a Elliptical carrier with ξ = 1.00 (RCP) @ λ = 14.5033 nm (T = 48.3777 as, ω = 3.1416 Ha = 85.4871 eV)
                                  – and a Gaussian envelope of duration 170.8811 as (intensity FWHM; ±2.00σ)
                                  – Uₚ = 0.0507 Ha = 1.3785 eV => α = 0.1433 Bohr = 7.5826 pm"""
        end
        withenv("UNITFUL_FANCY_EXPONENTS" => false) do
            @test string(F2) == """
                                Transversely polarized field with
                                  - I₀ = 2.0000e+00 au = 7.0188904e16 W cm^-2 =>
                                    - E₀ = 1.4142e+00 au = 727.2178 GV m^-1
                                    - A₀ = 0.4502 au
                                  – a Elliptical carrier with ξ = 1.00 (RCP) @ λ = 14.5033 nm (T = 48.3777 as, ω = 3.1416 Ha = 85.4871 eV)
                                  – and a Gaussian envelope of duration 170.8811 as (intensity FWHM; ±2.00σ)
                                  – Uₚ = 0.0507 Ha = 1.3785 eV => α = 0.1433 Bohr = 7.5826 pm"""
        end
    end

    @testset "Constant fields" begin
        @field(F) do
            I₀ = 4.0
            tmax = 4.0
            kind = :constant
        end

        withenv("UNITFUL_FANCY_EXPONENTS" => true) do
            @test string(F) == """
                               Constant field of
                                 - 4.0000 jiffies = 96.7554 as duration, and
                                 - E₀ = 2.0000e+00 au = 1.0284 TV m⁻¹"""
        end
        withenv("UNITFUL_FANCY_EXPONENTS" => false) do
            @test string(F) == """
                               Constant field of
                                 - 4.0000 jiffies = 96.7554 as duration, and
                                 - E₀ = 2.0000e+00 au = 1.0284 TV m^-1"""
        end

        @test F isa ElectricFields.ConstantField

        @test field_amplitude(F, -1) == 0
        @test field_amplitude(F, 2) ≈ 2 rtol=1e-7
        @test field_amplitude(F, 5) == 0

        @test intensity(F, -1) == 0
        @test intensity(F, 2) ≈ 4 rtol=1e-7
        @test intensity(F, 5) == 0

        @test field_amplitude(F, 0, 4.0) ≈ 8.0 rtol=1e-7

        @test intensity(F) ≈ 4.0 rtol=1e-7
        @test amplitude(F) ≈ 2.0 rtol=1e-7

        @test polarization(F) == LinearPolarization()

        # @test duration(F) == 4
        @test span(F) == 0..4

        @test isone(period(F))
        @test isone(max_frequency(F))
        @test photon_energy(F) == 2π

        @test dimensions(F) == 1

        @field(Fzero) do
            I₀ = 4.0
            tmax = 0.0
            kind = :constant
        end
        @test timeaxis(Fzero) == range(0, stop=0, length=0)
    end

    @testset "Ramps" begin
        @field(F) do
            I₀ = 4.0
            tmax = 4.0
            kind = :linear_ramp
        end

       @field(F2) do
            I₀ = 4.0
            tmax = 4.0
            kind = :linear_ramp
            ramp = :down
        end

        withenv("UNITFUL_FANCY_EXPONENTS" => true) do
            @test string(F) == """
                               Linear up-ramp of
                                 - 4.0000 jiffies = 96.7554 as duration, and
                                 - E₀ = 2.0000e+00 au = 1.0284 TV m⁻¹"""
            @test string(F2) == """
                                Linear down-ramp of
                                  - 4.0000 jiffies = 96.7554 as duration, and
                                  - E₀ = 2.0000e+00 au = 1.0284 TV m⁻¹"""
        end
        withenv("UNITFUL_FANCY_EXPONENTS" => false) do
            @test string(F) == """
                               Linear up-ramp of
                                 - 4.0000 jiffies = 96.7554 as duration, and
                                 - E₀ = 2.0000e+00 au = 1.0284 TV m^-1"""

            @test string(F2) == """
                                Linear down-ramp of
                                  - 4.0000 jiffies = 96.7554 as duration, and
                                  - E₀ = 2.0000e+00 au = 1.0284 TV m^-1"""
        end

        @test F isa ElectricFields.Ramp

        @test field_amplitude(F, -1) == 0
        @test field_amplitude(F, 0) ≈ 0 rtol=1e-7
        @test field_amplitude(F, 2) ≈ 1 rtol=1e-7
        @test field_amplitude(F, 4-2eps()) ≈ 2 rtol=1e-7
        @test field_amplitude(F, 5) == 0

        @test field_amplitude(F2, -1) == 0
        @test field_amplitude(F2, 0) ≈ 2 rtol=1e-7
        @test field_amplitude(F2, 2) ≈ 1 rtol=1e-7
        @test field_amplitude(F2, 4-2eps()) ≈ 0 atol=1e-10
        @test field_amplitude(F2, 5) == 0

        @test intensity(F, -1) == 0
        @test intensity(F, 0) ≈ 0 rtol=1e-7
        @test intensity(F, 2) ≈ 1 rtol=1e-7
        @test intensity(F, 4-2eps()) ≈ 4 rtol=1e-7
        @test intensity(F, 5) == 0

        @test intensity(F2, -1) == 0
        @test intensity(F2, 0) ≈ 4 rtol=1e-7
        @test intensity(F2, 2) ≈ 1 rtol=1e-7
        @test intensity(F2, 4-2eps()) ≈ 0 atol=1e-10
        @test intensity(F2, 5) == 0

        @test field_amplitude(F, 0, 4.0) ≈ 4.0 rtol=1e-7

        @test polarization(F) == LinearPolarization()

        # @test duration(F) == 4
        @test span(F) == 0..4

        @test isone(period(F))
        @test isone(max_frequency(F))
        @test photon_energy(F) == 2π

        @test dimensions(F) == 1

        @testset "Ramp kind = $(kind)" for (kind,half) in ((:linear_ramp, 0.5),
                                                           (:parabolic_ramp, 0.75),
                                                           (:sin²_ramp, 0.5))
            @field(R) do
                E₀ = 1.0
                tmax = 4.0u"fs"
                kind = kind
            end

            @field(R2) do
                E₀ = 1.0
                tmax = 4.0u"fs"
                kind = kind
                ramp = :down
            end

            @field(C) do
                E₀ = 1.0
                tmax = 4.0u"fs"
                kind = :constant
            end

            F = R + delay(C, duration(R))

            s = span(F)

            @test iszero(field_amplitude(F, s.left-1))
            @test iszero(field_amplitude(F, s.right+1))

            @test field_amplitude(R, austrip(2u"fs")) ≈ half
            @test field_amplitude(R2, austrip(2u"fs")) ≈ half

            @test vector_potential(F, s.right+1) ≈ vector_potential(F, s.right) rtol=1e-7

            t = range(0, stop=1.1span(F).right, length=1000)

            Fv = field_amplitude(F, t)
            Av = vector_potential(F, t)
            Fv2 = -ElectricFields.complex_derivative.(Ref(Base.Fix1(vector_potential, F)), t)

            @test Fv ≈ Fv2 rtol=1e-14
        end
    end
end


text:
"""
    add_reciprocal_links!(g, p, w)

Add reciprocal links to a weighted digraph, with probability `p` and weights sampled from `w`.

If `w` is a number, then it is assigned as the weight for every created link.
"""
function add_reciprocal_links!(g::MetaDiGraph, p::Real, w::Distribution = Uniform(0.0,1.0))
    for e in edges(g)
        (u, v) = src(e), dst(e)
        if !has_edge(g, v, u) && rand() < p
            add_edge!(g, v, u)
            set_prop!(g, v, u, :weight, rand(w))
        end
    end
    return g
end

function add_reciprocal_links!(g::MetaDiGraph, p::Real, w::Real)
    for e in edges(g)
        (u,v) = src(e), dst(e)
        if !has_edge(g, v, u) && rand() < p
            add_edge!(g, v, u)
            set_prop!(g, v, u, :weight, w)
        end
    end
    return g
end

"""
    add_reciprocal_links(g, p, w)

Add reciprocal links to a weighted digraph, with probability `p` and weights sampled from `w`.

This function creates a copy of `g`. See [add_reciprocal_links!](@ref) for the in-place version.
"""
add_reciprocal_links(g, p, w=Uniform(0.0,1.0)) = add_reciprocal_links!(copy(g), p, w)

add_reciprocal_links!(p::Real, w=Uniform(0.0,1.0)) = g -> add_reciprocal_links!(g, p, w)
add_reciprocal_links(p::Real, w=Uniform(0.0,1.0)) = g -> add_reciprocal_links(g, p, w)

ext: jl


text:
module Models

export Model, critical_point, local_boltzmann_weight, mpo, localham, checklocalham

struct Model{model} end

Model(s::AbstractString) = Model{Symbol(s)}()

# For notation:
# Model"tfim" == Model{:tfim}
# Model"heisenberg" == Model{:heisenberg}
macro Model_str(s)
  return :(Model{$(Expr(:quote, Symbol(s)))})
end

include("ising_classical_2d.jl")
include("hamiltonians.jl")

end


text:
""
function constraint_power_balance_shunt_dispatch(pm::AbstractACRModel, n::Int, i::Int, bus_arcs, bus_arcs_dc, bus_arcs_sw, bus_gens, bus_storage, bus_shunts_var, bus_pd, bus_qd, bus_gs_const, bus_bs_const)
    vi = var(pm, n, :vi, i)
    vr = var(pm, n, :vr, i)
    p    = get(var(pm, n),    :p, Dict()); PowerModels._check_var_keys(p, bus_arcs, "active power", "branch")
    q    = get(var(pm, n),    :q, Dict()); PowerModels._check_var_keys(q, bus_arcs, "reactive power", "branch")
    pg   = get(var(pm, n),   :pg, Dict()); PowerModels._check_var_keys(pg, bus_gens, "active power", "generator")
    qg   = get(var(pm, n),   :qg, Dict()); PowerModels._check_var_keys(qg, bus_gens, "reactive power", "generator")
    ps   = get(var(pm, n),   :ps, Dict()); PowerModels._check_var_keys(ps, bus_storage, "active power", "storage")
    qs   = get(var(pm, n),   :qs, Dict()); PowerModels._check_var_keys(qs, bus_storage, "reactive power", "storage")
    psw  = get(var(pm, n),  :psw, Dict()); PowerModels._check_var_keys(psw, bus_arcs_sw, "active power", "switch")
    qsw  = get(var(pm, n),  :qsw, Dict()); PowerModels._check_var_keys(qsw, bus_arcs_sw, "reactive power", "switch")
    p_dc = get(var(pm, n), :p_dc, Dict()); PowerModels._check_var_keys(p_dc, bus_arcs_dc, "active power", "dcline")
    q_dc = get(var(pm, n), :q_dc, Dict()); PowerModels._check_var_keys(q_dc, bus_arcs_dc, "reactive power", "dcline")

    bs = get(var(pm, n), :bs, Dict()); PowerModels._check_var_keys(bs, bus_shunts_var, "reactive power", "shunt")

    # possibly can save 2x in function eval, but no the dominant runtime at this moment
    #vm_sqr = @variable(pm.model, start=1.0, base_name="$(0)_vm_sqr_$(i)")

    #JuMP.@constraint(pm.model, vm_sqr == vi^2 + vr^2)
    #cstr_p = JuMP.@constraint(pm.model, 0 == - sum(p[a] for a in bus_arcs) + sum(pg[g] for g in bus_gens) - sum(pd for pd in values(bus_pd)) - sum(gs for gs in values(bus_gs_const))*vm_sqr)
    #cstr_q = JuMP.@constraint(pm.model, 0 == - sum(q[a] for a in bus_arcs) + sum(qg[g] for g in bus_gens) - sum(qd for qd in values(bus_qd)) + sum(bs for bs in values(bus_bs_const))*vm_sqr + sum(bs[s]*vm_sqr for s in bus_shunts_var))

    cstr_p = JuMP.@constraint(pm.model, 0 == - sum(p[a] for a in bus_arcs) + sum(pg[g] for g in bus_gens) - sum(pd for pd in values(bus_pd)) - sum(gs for gs in values(bus_gs_const))*(vi^2 + vr^2))
    cstr_q = JuMP.@NLconstraint(pm.model, 0 == - sum(q[a] for a in bus_arcs) + sum(qg[g] for g in bus_gens) - sum(qd for qd in values(bus_qd)) + sum(bs for bs in values(bus_bs_const))*(vi^2 + vr^2) + sum(bs[s]*(vi^2 + vr^2) for s in bus_shunts_var))

    if report_duals(pm)
        sol(pm, n, :bus, i)[:lam_kcl_r] = cstr_p
        sol(pm, n, :bus, i)[:lam_kcl_i] = cstr_q
    end
end

function constraint_ohms_yt_from_goc(pm::AbstractACRModel, n::Int, f_bus, t_bus, f_idx, t_idx, g, b, g_fr, b_fr, tr, ti, tm)
    p_fr = var(pm, n, :p, f_idx)
    q_fr = var(pm, n, :q, f_idx)
    vr_fr = var(pm, n, :vr, f_bus)
    vr_to = var(pm, n, :vr, t_bus)
    vi_fr = var(pm, n, :vi, f_bus)
    vi_to = var(pm, n, :vi, t_bus)

    JuMP.@constraint(pm.model, p_fr ==  (g/tm^2+g_fr)*(vr_fr^2 + vi_fr^2) + (-g*tr+b*ti)/tm^2*(vr_fr*vr_to + vi_fr*vi_to) + (-b*tr-g*ti)/tm^2*(vi_fr*vr_to - vr_fr*vi_to) )
    JuMP.@constraint(pm.model, q_fr == -(b/tm^2+b_fr)*(vr_fr^2 + vi_fr^2) - (-b*tr-g*ti)/tm^2*(vr_fr*vr_to + vi_fr*vi_to) + (-g*tr+b*ti)/tm^2*(vi_fr*vr_to - vr_fr*vi_to) )
end

""
function expression_branch_power_yt_from_goc(pm::AbstractACRModel, n::Int, f_bus, t_bus, f_idx, t_idx, g, b, g_fr, b_fr, tr, ti, tm)
    vr_fr = var(pm, n, :vr, f_bus)
    vr_to = var(pm, n, :vr, t_bus)
    vi_fr = var(pm, n, :vi, f_bus)
    vi_to = var(pm, n, :vi, t_bus)

    var(pm, n, :p)[f_idx] =  (g/tm^2+g_fr)*(vr_fr^2 + vi_fr^2) + (-g*tr+b*ti)/tm^2*(vr_fr*vr_to + vi_fr*vi_to) + (-b*tr-g*ti)/tm^2*(vi_fr*vr_to - vr_fr*vi_to)
    var(pm, n, :q)[f_idx] = -(b/tm^2+b_fr)*(vr_fr^2 + vi_fr^2) - (-b*tr-g*ti)/tm^2*(vr_fr*vr_to + vi_fr*vi_to) + (-g*tr+b*ti)/tm^2*(vi_fr*vr_to - vr_fr*vi_to)
end

""
function expression_branch_power_yt_from(pm::AbstractACRModel, n::Int, f_bus, t_bus, f_idx, t_idx, g, b, g_fr, b_fr, tr, ti, tm)
    vr_fr = var(pm, n, :vr, f_bus)
    vr_to = var(pm, n, :vr, t_bus)
    vi_fr = var(pm, n, :vi, f_bus)
    vi_to = var(pm, n, :vi, t_bus)

    var(pm, n, :p)[f_idx] =  (g+g_fr)/tm^2*(vr_fr^2 + vi_fr^2) + (-g*tr+b*ti)/tm^2*(vr_fr*vr_to + vi_fr*vi_to) + (-b*tr-g*ti)/tm^2*(vi_fr*vr_to - vr_fr*vi_to)
    var(pm, n, :q)[f_idx] = -(b+b_fr)/tm^2*(vr_fr^2 + vi_fr^2) - (-b*tr-g*ti)/tm^2*(vr_fr*vr_to + vi_fr*vi_to) + (-g*tr+b*ti)/tm^2*(vi_fr*vr_to - vr_fr*vi_to)
end

""
function expression_branch_power_yt_to(pm::AbstractACRModel, n::Int, f_bus, t_bus, f_idx, t_idx, g, b, g_to, b_to, tr, ti, tm)
    vr_fr = var(pm, n, :vr, f_bus)
    vr_to = var(pm, n, :vr, t_bus)
    vi_fr = var(pm, n, :vi, f_bus)
    vi_to = var(pm, n, :vi, t_bus)

    var(pm, n, :p)[t_idx] =  (g+g_to)*(vr_to^2 + vi_to^2) + (-g*tr-b*ti)/tm^2*(vr_fr*vr_to + vi_fr*vi_to) + (-b*tr+g*ti)/tm^2*(-(vi_fr*vr_to - vr_fr*vi_to))
    var(pm, n, :q)[t_idx] = -(b+b_to)*(vr_to^2 + vi_to^2) - (-b*tr+g*ti)/tm^2*(vr_fr*vr_to + vi_fr*vi_to) + (-g*tr-b*ti)/tm^2*(-(vi_fr*vr_to - vr_fr*vi_to))
end


text:
using CUDA, Test

function kernel_vadd(a, b, c)
    i = (blockIdx().x-1) * blockDim().x + threadIdx().x
    j = (blockIdx().y-1) * blockDim().y + threadIdx().y
    c[i,j] = a[i,j] + b[i,j]
    return nothing
end

@testset "Basic GPU addition" begin

    res = 1024

    # CUDAdrv functionality: generate and upload data
    a = round.(rand(Float32, (1024, 1024)) * 100)
    b = round.(rand(Float32, (1024, 1024)) * 100)
    d_a = CuArray(a)
    d_b = CuArray(b)
    d_c = similar(d_a)  # output array

    # run the kernel and fetch results
    # syntax: @cuda [kwargs...] kernel(args...)
    @cuda threads = (128, 1, 1) blocks = (div(res,128),res,1) kernel_vadd(d_a, d_b, d_c)

    # CUDAdrv functionality: download data
    # this synchronizes the device
    c = Array(d_c)
    a = Array(d_a)
    b = Array(d_b)

    @test a+b ≈ c
end

#GC.gc()

max_stars_repo_name: leios/simuleios


text:
module Examples

export greet, square, factorial

greet() = print("Hello World!")

function square(x)
    return x * x
end

function factorial(n)
    output = one(n)

    for number in 1:n
        output *= number
    end

    output  # The `return` is optional.
end

end # module

max_stars_repo_name: groundf/course-julia

max_issues_repo_name: groundf/course-julia

max_forks_repo_name: groundf/course-julia


text:
macro runtest(dir, prefix)
  quote
    fname = $dir * $prefix * ".jl"
    print("\n>>> Testing $fname\n\n")
    include(fname)
  end
end


text:
"""
    A container for information on aerosol size distribution
    and chemical properties.

    The size distribution is a sum of lognormal internally mixed modes.
    The chemical composition can be expressed using kappa parameter
    or hygroscopicity parameter B.
"""
module AerosolModel

export Mode_B
export Mode_κ

export AbstractAerosolDistribution

"""
    AbstractAerosolDistribution

The top-level super-type for all aerosol distribution types.
"""
abstract type AbstractAerosolDistribution{T} end

"""
    Mode_B

Represents the sizes and chemical composition
of aerosol particles in one size distribution mode.
The mode is assumed to be made up of internally mixed components
and follow a lognormal size distribution.
The chemical composition of aerosol particles in this mode
is described using the parameters from Abdul-Razzak and Ghan 2000.
"""
struct Mode_B{T, FT}
    "geometric mean dry radius"
    r_dry::FT
    "geometric standard deviation"
    stdev::FT
    "total number concentration"
    N::FT
    "tuple of mass mixing ratios for all components in this mode"
    mass_mix_ratio::T
    "tuple of mass fractions of soluble material for all components in this mode"
    soluble_mass_frac::T
    "tuple of osmotic coefficients for all components in this mode"
    osmotic_coeff::T
    "tuple of molar masses for all components in this mode"
    molar_mass::T
    "tuple of number of ions the salt dissociates into for all components in this mode"
    dissoc::T
    "tuple of aerosol densities for all components in this mode"
    aerosol_density::T
    "number of components in the mode"
    n_components::Int64
end

"""
    Mode_κ

Represents the sizes and chemical composition
of aerosol particles in one size distribution mode.
The mode is assumed to be made up of internally mixed components
and follow a lognormal size distribution.
The chemical composition of aerosol particles in this mode
is described using the parameters from Petters and Kreidenweis 2007.
"""
struct Mode_κ{T, FT}
    "geometric mean dry radius"
    r_dry::FT
    "geometric standard deviation"
    stdev::FT
    "total number concentration"
    N::FT
    "tuple of volume mixing ratios for all components in this mode"
    vol_mix_ratio::T
    "tuple of mass mixing ratios for all components in this mode"
    mass_mix_ratio::T
    "tuple of molar masses for all components in this mode"
    molar_mass::T
    "tuple of kappa-kohler values for all components in this mode"
    kappa::T
    "number of components in the mode"
    n_components::Int64
end

"""
    AerosolDistribution

Represents the aerosol size distribution as a tuple with different modes.
All modes have to either be of type Mode_B (Abdul-Razzak and Ghan 2000)
or of type Mode_κ (Petters and Kreidenweis 2007).

# Constructors

    AerosolDistribution(Modes::T)
"""
struct AerosolDistribution{T} <: AbstractAerosolDistribution{T}

    "tuple with all aerosol size distribution modes"
    Modes::T

    function AerosolDistribution(Modes::NTuple{N, T}) where {N, T}
        return new{typeof(Modes)}(Modes)
    end

end

end


text:
using Empirikos
using Test
using Random
using Hypatia

Z = BinomialSample(10)

betabin = marginalize(Z, Beta(2,1))

Zs = BinomialSample.(rand(MersenneTwister(1), betabin, 100), 10);
unique_Zs = sort(unique(Zs))
Zs_summary = summarize(Zs)

cvx_class = DiscretePriorClass(0.01:0.01:0.99)
npmle = NPMLE(cvx_class, Hypatia.Optimizer)
npmle_fit = fit(NPMLE(cvx_class, Hypatia.Optimizer), Zs)
npmle_fit_summary = fit(NPMLE(cvx_class, Hypatia.Optimizer), Zs_summary)

mosek_loglikelihood = -225.94642846739816
# figure out why tolerance has changed
@test loglikelihood(Zs, npmle_fit_summary.prior) ≈ mosek_loglikelihood atol = 0.1# 1e-3
@test loglikelihood(Zs, npmle_fit_summary.prior) ≈ loglikelihood(Zs, npmle_fit.prior) atol = 1e-3
@test loglikelihood(Zs, npmle_fit_summary.prior) ≈ loglikelihood(Zs_summary, npmle_fit_summary.prior) atol = 1e-3
@test loglikelihood(Zs, npmle_fit.prior) ≈ loglikelihood(Zs_summary, npmle_fit.prior) atol = 1e-3

@test probs(npmle_fit_summary.prior) ≈ probs(npmle_fit.prior) atol=1e-3

estimated_marginal_probs = pdf.(npmle_fit, unique_Zs)
true_marginal_probs = pdf.(betabin, response.(unique_Zs))

#plot(support(npmle_fit.prior), npmle_fit.prior, seriestype=:sticks)
#plot!(Beta(2,1))
#plot(0:10, true_marginal_probs, seriestype=:sticks)
#plot!(0:10, estimated_marginal_probs,  seriestype=:scatter)

marginalize.(unique_Zs, npmle_fit)

unique_binomial_samples = BinomialSample.(0:12, 12)
unique_binomial_counts = [3; 24; 104; 286; 670; 1033; 1343; 1112; 829; 478; 181; 45; 7]

dict_summary = Dict(unique_binomial_samples .=> unique_binomial_counts)

mult_binom = Empirikos.MultinomialSummary(dict_summary)
sort(collect(keys(mult_binom.store)))

beta_mle = fit(ParametricMLE(Beta()), mult_binom)

α_wiki = 34.1350
β_wiki = 31.6085
likelihood_wiki = -12492.9

@test beta_mle.α ≈ α_wiki atol = 0.05
@test beta_mle.β ≈ β_wiki atol = 0.05
@test loglikelihood(mult_binom, beta_mle) ≈ likelihood_wiki atol = 0.05


text:
@testset "647.palindromic-substrings.jl" begin
    @test count_substrings("abc") == 3
    @test count_substrings("aaa") == 6
    @test count_substrings("a" ^ 10) == sum(1:10)
    @test count_substrings("abcdefghijklmnnmlkjihgfedcba") == 42
end


text:
include("common/common.jl")
include("structured/structured.jl")


text:
const max_stimuluss=2.1
const epsilonn=1e-7
####################### c4 ####################

function c4_const(c4,c4_p)
    for i in 1:length(c4)
        c4[i]=0.1*c4_p[1]
    end
end

####################### c6 ####################

function c6_const(c6,c6_p)
    for i in 1:length(c6)
        c6[i]=c6_p[1]
    end
end

function c6_zeros(c6)
    for i in 1:length(c6)
        c6[i]=0.
    end
end

####################### c2 ####################

function c2_const(c2,c2_p)
    for i in 1:length(c2)
        c2[i]=c2_p[1]
    end
end

function c2_urgency(c2,c2_p)
    for i in 1:length(c2)
        c2[i]=c2_p[1]+0.1*c2_p[2]*(i-1)
    end
end

################# bias decision ##########

function bias_d_const(bias_d,bias_d_p)
    for i in 1:length(bias_d)
        bias_d[i]=bias_d_p[1]
    end
end

function bias_d_zeros(bias_d)
    for i in 1:length(bias_d)
        bias_d[i]=0.
    end
end

################# history bias decision ##########

function history_bias_d_const(hbias_d,bias_d_p)
    for i in 1:length(hbias_d)
        hbias_d[i]=bias_d_p[1]
    end
end

function history_bias_d_zeros(hbias_d)
    for i in 1:length(hbias_d)
        hbias_d[i]=0.
    end
end

##############################################3

function c2_kiani(c2,c2_p)
    for i in 1:length(c2)
        c2[i]=c2_p[i]
    end
end

###### initial condition#############################
function initial_condition(x0)
    #### initial condition
    #### x0_p[1] is the mean ###
    binN = length(x0)
    x0[Int( (binN-1)/2 )]=1
end

function initial_condition_bias(x0,bin_centers,x0_p)
    binN = length(x0)
    if x0_p[1]<bin_centers[1]
        x0[1]=1
    elseif  (x0_p[1]>bin_centers[end])
            x0[end]=1

    else
        if (x0_p[1] > bin_centers[1] && x0_p[1] < bin_centers[2])
            lp = 1; hp = 2;
        elseif (x0_p[1] > bin_centers[end-1] && x0_p[1] < bin_centers[end])
            lp = binN-1; hp = binN;
        else
            lp = floor(Int,((x0_p[1]-bin_centers[2])/dx)) + 2#find(bin_centers .<= sbins[k])[end]
            hp = ceil(Int,((x0_p[1]-bin_centers[2])/dx)) + 2#lp+1#Int(ceil((sbins[k]-bin_centers[2])/dx) + 1);
        end
        #println(k," ",lp," ",hp," ",sbin)
        if lp == hp
            x0[lp] = x0[lp] + 1
        else
            x0[hp] = x0[hp] + 1*(x0_p[1] - bin_centers[lp])/(bin_centers[hp] - bin_centers[lp])
            x0[lp] = x0[lp] + 1*(bin_centers[hp] - x0_p[1])/(bin_centers[hp] - bin_centers[lp])
        end
    end

end

function initial_condition_bias_hbias(x0,bin_centers,bias_total)

    #Overall bias x0_p[1]
    #hbias_x0 parameter x0_p[2]
    #historybias from data x0_p[3]

    Ntrials,binN = size(x0)

    for itrial in 1:Ntrials

        if bias_total[itrial]<bin_centers[1]
            x0[itrial,1]=1
        elseif  (bias_total[itrial]>bin_centers[end])
                x0[itrial,end]=1

        else
            if (bias_total[itrial] > bin_centers[1] && bias_total[itrial] < bin_centers[2])
                lp = 1; hp = 2;
            elseif (bias_total[itrial] > bin_centers[end-1] && bias_total[itrial] < bin_centers[end])
                lp = binN-1; hp = binN;
            else
                lp = floor(Int,((bias_total[itrial]-bin_centers[2])/dx)) + 2#find(bin_centers .<= sbins[k])[end]
                hp = ceil(Int,((bias_total[itrial]-bin_centers[2])/dx)) + 2#lp+1#Int(ceil((sbins[k]-bin_centers[2])/dx) + 1);
            end
            #println(k," ",lp," ",hp," ",sbin)
            if lp == hp
                x0[itrial,lp] = x0[lp] + 1
            else
                x0[itrial,hp] = x0[itrial,hp] + 1*(bias_total[itrial] - bin_centers[lp])/(bin_centers[hp] - bin_centers[lp])
                x0[itrial,lp] = x0[itrial,lp] + 1*(bin_centers[hp] - bias_total[itrial])/(bin_centers[hp] - bin_centers[lp])
            end
        end
    end
end

function initial_condition_bias_variance(x0,bin_centers,x0_p)
    #### we add variability and bias in the initial condition
    #### x0_p[1] is the mean, x0_p[2] is the std

    binN = length(bin_centers)
    sigma0=x0_p[2]
    n_sbins = 11#max(70, ceil(10*x0_p[2]/dx))
    swidth = 5.0*sigma0
    println("hola")
    println(swidth)
    println(n_sbins)
    sbinsize = swidth/n_sbins #
    base_sbins    = collect(-swidth:sbinsize:swidth)
    ps= exp.(-base_sbins.^2/(2*(x0_p[2]^2)))
    ps= ps/sum(ps);
    sbin_length = length(base_sbins)
    for k in 1:sbin_length
        sbin = (k-1)*sbinsize + x0_p[1] - swidth
        #println(sbin)
        if sbin <= bin_centers[1] #(bin_centers[1] + bin_centers[2])/2
            x0[1] = x0[1] + ps[k]
        elseif bin_centers[end] <= sbin#(bin_centers[end]+bin_centers[end-1])/2 <= sbins[k]
            x0[end] = x0[end] + ps[k]
        else
            if (sbin > bin_centers[1] && sbin < bin_centers[2])
                lp = 1; hp = 2;
            elseif (sbin > bin_centers[end-1] && sbin < bin_centers[end])
                lp = binN-1; hp = binN;
            else
                lp = floor(Int,((sbin-bin_centers[2])/dx)) + 2#find(bin_centers .<= sbins[k])[end]
                hp = ceil(Int,((sbin-bin_centers[2])/dx)) + 2#lp+1#Int(ceil((sbins[k]-bin_centers[2])/dx) + 1);
            end
            #println(k," ",lp," ",hp," ",sbin)
            if lp == hp
                x0[lp] = x0[lp] + ps[k]
            else
                x0[hp] = x0[hp] + ps[k]*(sbin - bin_centers[lp])/(bin_centers[hp] - bin_centers[lp])
                x0[lp] = x0[lp] + ps[k]*(bin_centers[hp] - sbin)/(bin_centers[hp] - bin_centers[lp])
            end
        end
    end

    return x0
end

#######################stimulus transformation####################
function st_linear(s,s_trans,st_p)
    for iframe in 1:length(s)
        s_trans[iframe]=s[iframe]*(st_p[1])
    end
end

function st_linear_adaptation(s,s_trans,st_p)
    for iframe in 1:length(s)
        s_trans[iframe]=s[iframe]*(st_p[1]-0.1*st_p[2]*(iframe-1) )
    end
end

function st_linear_adaptation_plot(s,s_trans,st_p,iframe)
    for i in 1:length(s)
        s_trans[i]=s[i]*(st_p[1]-0.1*st_p[2]*(iframe-1) )
    end
end

function st_linear_adaptation_exp(s,s_trans,st_p)
    for iframe in 1:length(s)
        s_trans[iframe]=s[iframe]*(st_p[1]*exp(st_p[2]*(iframe-1)))
    end
end

function st_linear_adaptation_exp_bias(s,s_trans,st_p)
    for iframe in 1:length(s)
        s_trans[iframe]=s[iframe]*(st_p[1]*exp(st_p[2]*(iframe-1)) -st_p[3])
    end
end

function st_linear_adaptation_exp_plot(s,s_trans,st_p,iframe)
    for i in 1:length(s)
        s_trans[i]=s[i]*(st_p[1]*exp(st_p[2]*(iframe-1)) )
    end
end

function st_logistic_frame(s,st_p,k)
    beta= sign(st_p[2])*maximum([epsilonn,abs(st_p[2]) ]) +epsilonn/10.0
    return (st_p[1]/k)*(1.0/(1+exp(-( beta*s  ) ) )-0.5 )

end

function st_logistic(stim,s_trans,st_p)
    beta= sign(st_p[2])*maximum([epsilonn,abs(st_p[2]) ]) +epsilonn/10.0
    println(beta)
    k=1.0/(1+exp(-beta*max_stimuluss) ) -0.5

    Nframes=length(stim)
    for iframe in 1:Nframes
        #s_trans[iframe]=(st_p[1]/k)*(1.0/(1+exp(-( beta*stim[iframe]  ) ) )-0.5 )
        s_trans[iframe]=(1.0/(1+exp(-( beta*stim[iframe]  ) ) )-0.5 )

    end
end

function st_logistic_bias(stim,s_trans,st_p)
    beta= sign(st_p[2])*maximum([epsilonn,abs(st_p[2]) ]) +epsilonn/10.0
    println(beta)
    k=1.0/(1+exp(-beta*max_stimuluss) ) -0.5

    Nframes=length(stim)
    for iframe in 1:Nframes
        if stim[iframe]>0
            s_trans[iframe]=st_logistic_frame(stim[iframe],[st_p[1],st_p[2]],k)
        else
            s_trans[iframe]=st_logistic_frame(stim[iframe],[st_p[1]-st_p[3],st_p[2]],k)
        end
    end
end

function st_logistic_adaptation(stim,s_trans,st_p)
    beta= sign(st_p[2])*maximum([epsilonn,abs(st_p[2]) ]) +epsilonn/10.0
    k=1.0/(1+exp(-beta*max_stimuluss) ) -0.5

    Nframes=length(stim)
    for iframe in 1:Nframes
        s_trans[iframe]=st_logistic_frame(stim[iframe],[st_p[1]-0.1*st_p[3]*(iframe-1),st_p[2]],k)
    end
end

function st_logistic_adaptation_plot(stim,s_trans,st_p,ifrmae)

    k=1.0/(1+exp(-beta*max_stimuluss) ) -0.5

    Nframes=length(stim)
    for i in 1:Nframes
        s_trans[i]=st_logistic_frame(stim[i],[st_p[1]-0.1*st_p[3]*(iframe-1),st_p[2]],k)
    end
end

function st_logistic_adaptation_bias(stim,s_trans,st_p)

    beta= sign(st_p[2])*maximum([epsilonn,abs(st_p[2]) ]) +epsilonn/10.0
    k=1.0/(1+exp(-beta*max_stimuluss) ) -0.5

    Nframes=length(stim)
    for iframe in 1:Nframes
        if stim[iframe]>0
            s_trans[iframe]=st_logistic_frame(stim[iframe],[st_p[1]-0.1*st_p[3]*(iframe-1),st_p[2]],k)
        else
            s_trans[iframe]=st_logistic_frame(stim[iframe],[st_p[1]-0.1*st_p[3]*(iframe-1)-st_p[4],st_p[2]],k)
        end
    end
end

function st_logistic_adaptation_bias_plot(stim,s_trans,st_p,iframe)

    beta= sign(st_p[2])*maximum([epsilonn,abs(st_p[2]) ]) +epsilonn/10.0
    k=1.0/(1+exp(-beta*max_stimuluss) ) -0.5

    Nframes=length(stim)
    for i in 1:Nframes
        if stim[iframe]>0
            s_trans[i]=st_logistic_frame(stim[i],[st_p[1]-0.1*st_p[3]*(iframe-1),st_p[2]],k)
        else
            s_trans[i]=st_logistic_frame(stim[i],[st_p[1]-0.1*st_p[3]*(iframe-1)-st_p[4],st_p[2]],k)
        end
    end
end

# function st_logistic(stim,s_trans,st_p)
#     Nframes=length(stim)
#     for iframe in 1:Nframes
#         s_trans[iframe]=2.0*(1/st_p[1])*(1.0/(1+exp(-( st_p[2]*st_p[1]*stim[iframe]  ) ) )-0.5 )
#     end
# end

#
# function st_logistic_adaptation(stim,s_trans,st_p)
#     Nframes=length(stim)
#     for iframe in 1:Nframes
#         s_trans[iframe]=2.0*(1/st_p[1])*(1.0/(1+exp(-( (st_p[2]-0.1*st_p[3]*(iframe-1))*st_p[1]*stim[iframe]  ) ) )-0.5 )
#         #s_trans[iframe]=2.0*(1.0/st_p[1])*(1.0/(1+exp(-( (1/(st_p[2]-0.1*st_p[3]*(iframe-1)) )*stim[iframe] ) ) )-0.5 )
#     end
# end
#
#
# function st_logistic_adaptation_plot(stim,s_trans,st_p,iframe)
#     Nframes=length(stim)
#     for i in 1:Nframes
#         s_trans[i]=2.0*(1/st_p[1])*(1.0/(1+exp(-( (st_p[2]-0.1*st_p[3]*(iframe-1))*st_p[1]*stim[i]  ) ) )-0.5 )
#     end
# end
#

#
# function st_logistic_bias(stim,s_trans,st_p)
#     Nframes=length(stim)
#     for iframe in 1:Nframes
#             if stim[iframe]-st_p[3] > 0
#                 s_trans[iframe]=2*st_p[1]*(1.0/(1+exp(-(st_p[2]*(stim[iframe]-st_p[3] ) ) ) ) -0.5 )
#             else
#                 #s_trans[iframe]=2*st_p[1]*(1./(1+exp(-(st_p[2]*(stim[iframe]-st_p[3] ) ) ) ) -0.5 )
#                 s_trans[iframe]=2*(st_p[1]-st_p[4])*(1.0/(1+exp(-( (st_p[2]-st_p[5])*(stim[iframe]-st_p[3] ) ) ) ) -0.5 )
#             end
#     end
# end

# function st_logistic_bias(stim,s_trans,st_p)
#     Nframes=length(stim)
#     for iframe in 1:Nframes
#             if stim[iframe] > 0
#                 s_trans[iframe]=2.0*(1/st_p[1])*(1.0/(1+exp(-( st_p[2]*st_p[1]*stim[iframe]  ) ) )-0.5 )
#
#                 #s_trans[iframe]=2*(1.0/st_p[1])*(1.0/(1+exp(-( (1./st_p[2])*(stim[iframe] ) ) ) ) -0.5 )
#             else
#                 #s_trans[iframe]=2*st_p[1]*(1./(1+exp(-(st_p[2]*(stim[iframe]-st_p[3] ) ) ) ) -0.5 )
#                 #s_trans[iframe]=2*( (1.0/st_p[1])-st_p[3])*(1.0/(1+exp(-( ( (1./st_p[2])-st_p[4])*(stim[iframe] ) ) ) ) -0.5 )
#                 s_trans[iframe]=2.0*( 1/st_p[1] -st_p[3] )*(1.0/(1+exp(-( (st_p[2]-st_p[4])*st_p[1]*stim[iframe]  ) ) )-0.5 )
#             end
#     end
# end
#
#
#
# function st_logistic_adaptation_bias(stim,s_trans,st_p)
#     Nframes=length(stim)
#     for iframe in 1:Nframes
#             if stim[iframe] > 0
#                 #s_trans[iframe]=2*(1.0/st_p[1])*(1.0/(1+exp(-( (1./(st_p[2]-0.1*st_p[3]*(iframe-1))    )*(stim[iframe] ) ) ) ) -0.5 )
#                 s_trans[iframe]=2.0*(1/st_p[1])*(1.0/(1+exp(-( (st_p[2]-0.1*st_p[3]*(iframe-1))*st_p[1]*stim[iframe] )  ) )-0.5 )
#
#             else
#                 #s_trans[iframe]=2*( (1.0/st_p[1])-st_p[4])*(1.0/(1+exp(-( ( (1./ (st_p[2]-0.1*st_p[3]*(iframe-1))  )-st_p[5])*(stim[iframe] ) ) ) ) -0.5 )
#                 s_trans[iframe]=2.0*(1/st_p[1]-st_p[4])*(1.0/(1+exp(-( ( (st_p[2]-0.1*st_p[3]*(iframe-1)) -st_p[5] )*st_p[1]*stim[iframe] )  ) )-0.5 )
#
#             end
#     end
# end
#
#
#
#
# function st_logistic_adaptation_bias_plot(stim,s_trans,st_p,iframe)
#     Nframes=length(stim)
#     for i in 1:Nframes
#             if stim[i] > 0
#                 #s_trans[iframe]=2*(1.0/st_p[1])*(1.0/(1+exp(-( (1./(st_p[2]-0.1*st_p[3]*(iframe-1))    )*(stim[iframe] ) ) ) ) -0.5 )
#                 s_trans[i]=2.0*(1/st_p[1])*(1.0/(1+exp(-( (st_p[2]-0.1*st_p[3]*(iframe-1))*st_p[1]*stim[i]   ) ) ) -0.5 )
#
#             else
#                 #s_trans[iframe]=2*( (1.0/st_p[1])-st_p[4])*(1.0/(1+exp(-( ( (1./ (st_p[2]-0.1*st_p[3]*(iframe-1))  )-st_p[5])*(stim[iframe] ) ) ) ) -0.5 )
#                 s_trans[i]=2.0*(1/st_p[1]-st_p[4])*(1.0/(1+exp(-( ( (st_p[2]-0.1*st_p[3]*(iframe-1)) -st_p[5] )*st_p[1]*stim[i] ) ) )-0.5 )
#
#             end
#     end
# end
#

##############sigma ########
function sigma_const(sigma,sigma_p)
    for iframe in 1:length(sigma)
        sigma[iframe]=sigma_p[1]
    end
    #return s_trans
end

function sigma_delay_kiani(sigma,sigma_p)
    for iframe in 1:length(sigma)
        if iframe%2==0
            sigma[iframe]=sigma_p[1]/sqrt(2.)
        else
            sigma[iframe]=sigma_p[1]
        end
    end
    #return s_trans
end

#######################dv  read out####################

function read_out_soft(read_out,bins,ro_p)
    for ibin in 1:length(bins)
        read_out[ibin]=1.0/(1+exp.(-10*ro_p[1]*(bins[ibin]-ro_p[2])))
    end

end

function read_out_perfect(read_out,binN)
    read_out[1:binN]=0
    read_out[binN+1]=0.5
    read_out[binN+2:end]=1
end

function potential_DW(param,x)
    return -param[0]-param[1]*x+pram[2]*x^3
end

function read_out_TW(ro,bin_centers,c2,c4,c6)
    delta=c4^2+4*c2*c6
    if  (c2<0) & (c4<0) & (c6>0) & (delta>0)
        b=sqrt( (-c4-sqrt(delta))/(2*c6)   )
        i=1
        accu=0.
        while (i<=length(ro)) & (bin_centers[i]<-b)
            ro[i]=0
            i=i+1
        end
        while (i<=length(ro)) & (bin_centers[i]<b)
            ro[i]=0.5
            i=i+1
        end
        while i<=length(ro)
            ro[i]=1
            i=i+1
        end
    elseif (c2<0) & (c4>0) & (c6>0)
        ro[:]=0.5

    else
        N=(length(bin_centers)-1)/2
        ro[:]=vcat(zeros(N),[0.5],ones(N))
    end
    return ro
end

function read_out_DW(ro,bin_centers,c2,c4)
    if  (c2<0) & (c4<0)

        b=sqrt(c2/c4)
        i=1
        accu=0.
        while  (i<=length(ro)) & (bin_centers[i]<-b)
            ro[i]=0
            i=i+1
        end
        while  (i<=length(ro)) & (bin_centers[i]<-b)
            ro[i]=0.5
            i=i+1
        end
        while i<=length(ro)
            ro[i]=1
            i=i+1
        end
    elseif (c2<0) & (c4>0)
        ro[:]=0.5

    else
        N=(length(bin_centers)-1)/2
        ro[:]=vcat(zeros(N),[0.5],ones(N))
    end
    return ro
end

max_stars_repo_name: emmafdz/FinalDegreeProject

max_issues_repo_name: emmafdz/FinalDegreeProject

max_forks_repo_name: emmafdz/FinalDegreeProject


text:
using SafeTestsets

@safetestset "Vulkan.jl" begin
    include("base.jl")
end

# Integration tests on unregistered packages.
# It is advised not to try to install them because
# they have also unregistered dependencies which makes
# it very painful to setup.

const DEV_PATH = joinpath(homedir(), ".julia", "dev")

using Pkg

if isdir(DEV_PATH)
    for package in ["Givre", "VulkanExamples"]
        if package in readdir(DEV_PATH)
            Pkg.activate(joinpath(DEV_PATH, package))
            Pkg.test()
            GC.gc()
        end
    end
end


text:
struct LSSymbolServerFailure <: Exception
    msg::AbstractString
end

function Base.showerror(io::IO, ex::LSSymbolServerFailure)
    print(io, ex.msg)
end

struct LSUriConversionFailure <: Exception
    msg::AbstractString
end

function Base.showerror(io::IO, ex::LSUriConversionFailure)
    print(io, ex.msg)
end

struct LSOffsetError <: Exception
    msg::AbstractString
end

function Base.showerror(io::IO, ex::LSOffsetError)
    print(io, ex.msg)
end

struct LSSyncMismatch <: Exception
    msg::AbstractString
end

function Base.showerror(io::IO, ex::LSSyncMismatch)
    print(io, ex.msg)
end

struct LSHoverError <: Exception
    msg::AbstractString
end

function Base.showerror(io::IO, ex::LSHoverError)
    print(io, ex.msg)
end

struct LSPositionToOffsetException <: Exception
    msg::AbstractString
end

function Base.showerror(io::IO, ex::LSPositionToOffsetException)
    print(io, ex.msg)
end

struct LSRelativePath <: Exception
    msg::AbstractString
end

function Base.showerror(io::IO, ex::LSRelativePath)
    print(io, ex.msg)
end

struct LSInfiniteLoop <: Exception
    msg::AbstractString
end

function Base.showerror(io::IO, ex::LSInfiniteLoop)
    print(io, ex.msg)
end

struct LSInvalidFile <: Exception
    msg::AbstractString
end

function Base.showerror(io::IO, ex::LSInvalidFile)
    print(io, ex.msg)
end


text:
function bubblesort!(arr::AbstractVector)
    for _ in 2:length(arr), j in 1:length(arr)-1
        if arr[j] > arr[j+1]
            arr[j], arr[j+1] = arr[j+1], arr[j]
        end
    end
    return arr
end

v = rand(-10:10, 10)
println("# unordered: $v\n -> ordered: ", bubblesort!(v))

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
using Test
using SPGBox

@testset "simple polynomial" begin

  f(x) = x[1]^2 + (x[2]-1.)^2

  function g!(g,x)
    g[1] = 2*x[1]
    g[2] = 2*(x[2]-1.)
  end

  x = [ 10. , 18. ]
  R = spgbox!(f,g!,x)
  @test R.f ≈ 0.
  @test R.x ≈ [0.,1.]

  x = [ 10. , 18. ]
  R = spgbox!(f,g!,x,lower=[2,-Inf])
  @test R.f ≈ 4.
  @test R.x ≈ [2.,1.]

  x = [ 10. , 18. ]
  R = spgbox!(f,g!,x,lower=[-Inf,2.])
  @test R.f ≈ 1.
  @test R.x ≈ [0.,2.]

  x = [ 10. , 18. ]
  R = spgbox!(f,g!,x,upper=[+Inf,-2])
  @test R.f ≈ 9.
  @test R.x ≈ [0.,-2.]

  x = [ 10. , 18. ]
  R = spgbox!(f,g!,x,upper=[-2,+Inf])
  @test R.f ≈ 4.
  @test R.x ≈ [-2.,1.]

  # just testing the interface
  x = [ 10. , 18. ]
  R = spgbox!(f,g!,x,lower=[-Inf,2.])
  @test R.f ≈ 1.
  @test R.x ≈ [0.,2.]

  x = [ 10. , 18. ]
  R = spgbox!(f,g!,x,upper=[-2,+Inf])
  @test R.f ≈ 4.
  @test R.x ≈ [-2.,1.]

  x = [ 10. , 18. ]
  R = spgbox!(f,g!,x,lower=[-Inf,-2.],upper=[+Inf,2])
  @test R.f ≈ 0.
  @test R.x ≈ [0.,1.]

  x = [ 10. , 18. ]
  R = spgbox!(f,g!,[-Inf,-2.],[+Inf,2],x)
  @test R.f ≈ 0.
  @test R.x ≈ [0.,1.]

  # Test the mutating call
  x = [ 10. , 18. ]
  R = spgbox(f,g!,[-Inf,-2.],[+Inf,2],x)
  @test x == [ 10. , 18. ]
  @test R.f ≈ 0.
  @test R.x ≈ [0.,1.]

  x = [ 10. , 18. ]
  R = spgbox(f,g!,x,lower=[-Inf,-2.],upper=[+Inf,2])
  @test x == [ 10. , 18. ]
  @test R.f ≈ 0.
  @test R.x ≈ [0.,1.]

end


text:
module LazyImports

export @shim_import, @lazy_import

using Base: Docs
using Pkg: TOML
using UUIDs

abstract type AbstractModuleShim end

struct ShimModule <: AbstractModuleShim
    uuid::UUID
    name::Symbol
end

struct LazyModule <: AbstractModuleShim
    uuid::UUID
    name::Symbol
end

_uuid(m::AbstractModuleShim) = getfield(m, :uuid)
Base.nameof(m::AbstractModuleShim) = getfield(m, :name)

Base.PkgId(m::AbstractModuleShim) = Base.PkgId(_uuid(m), String(nameof(m)))

getmodule(m::ShimModule) = Base.loaded_modules[Base.PkgId(m)]
getmodule(m::LazyModule) = Base.require(Base.PkgId(m))

Base.getproperty(m::AbstractModuleShim, name::Symbol) =
    getproperty(getmodule(m), name)

Docs.Binding(m::AbstractModuleShim, name::Symbol) =
    Docs.Binding(getmodule(m), name)

topmodule(m::Module) = parentmodule(m) == m ? m : topmodule(parentmodule(m))

function _uuidfor(__module__::Module, name::Symbol)
    root = topmodule(__module__)
    jlpath = pathof(root)
    if jlpath == nothing
        error("""
        Module $__module__ does not have associated file.
        """)
    end
    tomlpath = joinpath(dirname(dirname(jlpath)), "Project.toml")
    if !isfile(tomlpath)
        error("""
        `LazyImports` needs package `$(nameof(root))` to use `Project.toml`
        file.  `Project.toml` is not found at:
            $tomlpath
        """)
    end

    found = find_uuid_or(TOML.parsefile(tomlpath), String(name)) do
        error("""
Package `$name` is not listed in `[deps]` or `[extras]` of `Project.toml`
file for `$(nameof(root))` found at:
    $tomlpath
If you are developing `$(nameof(root))`, add `$name` to the dependency.
Otherwise, please report this to `$(nameof(root))`'s issue tracker.
""")
    end

    # Just to be extremely careful, editing Project.toml should
    # invalidate the compilation cache since the UUID may be changed
    # or removed:
    include_dependency(tomlpath)

    return UUID(found)
end

find_uuid_or(f, project::Dict, name::String) =
    get(get(project, "deps", Dict()), name) do
        get(get(project, "extras", Dict()), name) do
            f()
        end
    end

_uuid(id::UUID) = id
_uuid(id) = UUID(id)

create(factory, __module__, name) = factory(_uuidfor(__module__, name), name)

import_expr(factory, __module__, name::Symbol) =
    :(const $name = $create($factory, $__module__, $(QuoteNode(name))))

function import_expr(factory, ::Any, expr::Expr)
    @assert expr.head == :(=)
    name, uuid = expr.args
    :(const $name = $factory($_uuid($uuid), $(QuoteNode(name))))
end

"""
    @shim_import Module
    @shim_import Module = UUID

This acts like `@lazy_import` except that `Module.f` works only if
`Module` is already loaded.
"""
macro shim_import(expr)
    esc(import_expr(ShimModule, __module__, expr))
end

"""
    @lazy_import Module
    @lazy_import Module = UUID

Construct a module-like object for `Module`.  Accessing `Module.f`
imports `Module` and then get `f` from `Module`.  The `UUID` can be
omitted if `Module` is listed in `[deps]` or `[extras]` of
`Project.toml` file of the package using `@lazy_import`.
"""
macro lazy_import(expr)
    esc(import_expr(LazyModule, __module__, expr))
end

end # module


text:
# one-body non-bonded potentials
function Singlea(xi, qi, ti, mu, eta, kappa)
    u = 0.0;
    return u;
end

# one-body bonded potentials
function Singleb(xi, qi, ti, mu, eta, kappa)
    u = 0.0;
    return u;
end

function ezbl(r,zi,zj)
  pzbl = 0.23;
  a0 = 0.46850;
  c1 = 0.02817;
  c2 = 0.28022;
  c3 = 0.50986;
  c4 = 0.18175;
  d1 = 0.20162;
  d2 = 0.40290;
  d3 = 0.94229;
  d4 = 3.19980;
  qqr2e = 14.399645;
  qelectron = 1.0;

  ainv = (sympy.Pow(zi,pzbl) + sympy.Pow(zj,pzbl))/(a0);
  d1a = d1*ainv;
  d2a = d2*ainv;
  d3a = d3*ainv;
  d4a = d4*ainv;
  zze = zi*zj*qqr2e*qelectron*qelectron;

  #[d1a d2a d3a d4a] #-> 2.3090    4.6140   10.7912   36.6444

  e1 = exp(-d1a*r);
  e2 = exp(-d2a*r);
  e3 = exp(-d3a*r);
  e4 = exp(-d4a*r);

  rinv = 1/r;
  sum = c1*e1;
  sum = sum + c2*e2;
  sum = sum + c3*e3;
  sum = sum + c4*e4;

  sum_p = -c1*d1a*e1;
  sum_p = sum_p - c2*d2a*e2;
  sum_p = sum_p - c3*d3a*e3;
  sum_p = sum_p - c4*d4a*e4;

  sum_pp = c1*e1*d1a*d1a;
  sum_pp = sum_pp + c2*e2*d2a*d2a;
  sum_pp = sum_pp + c3*e3*d3a*d3a;
  sum_pp = sum_pp + c4*e4*d4a*d4a;

  e = zze*sum*rinv;
  ep = zze*(sum_p - sum*rinv)*rinv;
  epp = zze*(sum_pp - 2.0*sum_p*rinv + 2.0*sum*rinv*rinv)*rinv;

  return e, ep, epp
end

# two-body non-bonded potentials
function Paira(xij, qi, qj, ti, tj, mu, eta, kappa)

    r2 = xij[1]*xij[1] + xij[2]*xij[2] + xij[3]*xij[3];
    r = sqrt(r2);

    cut_inner = 4;
    cut_outer = 4.8;
    z1 = 74;
    z2 = 4;

    # (zi, zj) parameters depend on atom types
    zi = z1*(2 - ti[1]) + z2*(ti[1]-1); # ti = 1 => zi = z1, ti = 2 => zi = z2
    zj = z1*(2 - tj[1]) + z2*(tj[1]-1); # tj = 1 => zj = z1, tj = 2 => zj = z2

    # ZBL potential
    u = ezbl(r,zi,zj)[1];

    tc = cut_outer - cut_inner;
    fc, fcp, fcpp = ezbl(cut_outer, zi, zj);
    swa = (-3.0*fcp + tc*fcpp)/(tc*tc);
    swb = ( 2.0*fcp - tc*fcpp)/(tc*tc*tc);
    swc = -fc + (tc/2.0)*fcp - (tc*tc/12.0)*fcpp;

    #sw1 = swa;
    #sw2 = swb;
    sw3 = swa/3.0;
    sw4 = swb/4.0;
    sw5 = swc;
    #[sw3 sw4 sw5] # ->  [0.0456   -0.0343   -0.0163]

    t = r - cut_inner;
    sinner = sw5;
    souter = sw5 + t*t*t * (sw3 + sw4*t);

    # switching function
    a = 0.5 + 0.5*tanh(1e2*t); # a = 0 when t<0, and a = 1 when t>0
    s = a*souter + (1-a)*sinner;
    u = u + s;
    return u;
end

# two-body bonded potentials
function Pairb(xij, qi, qj, ti, tj, mu, eta, kappa)
    u = 0.0;
    return u;
end

# two-body bond order potentials (EAM-like potentials)
function Pairc(xij, qi, qj, ti, tj, rho, mu, eta, kappa)
    u = 0.0;
    return u;
end

# three-body non-bonded potentials
function Tripleta(xij, xik, qi, qj, qk, ti, tj, tk, mu, eta, kappa)

    u = 0.0;
    return u;
end

# three-body bonded potentials
function Tripletb(xij, xik, qi, qj, qk, ti, tj, tk, mu, eta, kappa)
    u = 0.0;
    return u;
end

# three-body bond porder potentials (Terssof-like potentials)
function Tripletc(xij, xik, qi, qj, qk, ti, tj, tk, rho, mu, eta, kappa)
    u = 0.0;
    return u;
end

# four-body non-bonded potentials
function Quadrupleta(xij, xik, xil, qi, qj, qk, ql, ti, tj, tk, tl, mu, eta, kappa)
    u = 0.0;
    return u;
end

# four-body bonded potentials
function Quadrupletb(xij, xik, xil, qi, qj, qk, ql, ti, tj, tk, tl, mu, eta, kappa)
    u = 0.0;
    return u;
end

max_stars_repo_name: cesmix-mit/MDP

max_issues_repo_name: cesmix-mit/MDP

max_forks_repo_name: cesmix-mit/MDP


text:
# This file is a part of Julia. License is MIT: http://julialang.org/license

module LibGit2

import Base: merge!, cat, ==

export with, GitRepo, GitConfig

const GITHUB_REGEX =
    r"^(?:git@|git://|https://(?:[\w\.\+\-]+@)?)github.com[:/](([^/].+)/(.+?))(?:\.git)?$"i

include("utils.jl")
include("consts.jl")
include("types.jl")
include("error.jl")
include("signature.jl")
include("oid.jl")
include("reference.jl")
include("commit.jl")
include("repository.jl")
include("config.jl")
include("walker.jl")
include("remote.jl")
include("strarray.jl")
include("index.jl")
include("merge.jl")
include("tag.jl")
include("blob.jl")
include("diff.jl")
include("rebase.jl")
include("status.jl")
include("tree.jl")
include("callbacks.jl")

using .Error

immutable State
    head::Oid
    index::Oid
    work::Oid
end

"""Return HEAD Oid as string"""
function head(pkg::AbstractString)
    with(GitRepo, pkg) do repo
        string(head_oid(repo))
    end
end

""" git update-index """
function need_update(repo::GitRepo)
    if !isbare(repo)
        # read updates index from filesystem
        read!(repo, true)
    end
end

""" Checks if commit is in repository """
function iscommit(id::AbstractString, repo::GitRepo)
    res = true
    try
        c = get(GitCommit, repo, id)
        if c === nothing
            res = false
        else
            finalize(c)
        end
    catch
        res = false
    end
    return res
end

""" git diff-index HEAD [-- <path>]"""
isdirty(repo::GitRepo, paths::AbstractString=""; cached::Bool=false) =
    isdiff(repo, Consts.HEAD_FILE, paths, cached=cached)

""" git diff-index <treeish> [-- <path>]"""
function isdiff(repo::GitRepo, treeish::AbstractString, paths::AbstractString=""; cached::Bool=false)
    tree_oid = revparseid(repo, "$treeish^{tree}")
    iszero(tree_oid) && return true
    result = false
    tree = get(GitTree, repo, tree_oid)
    try
        diff = diff_tree(repo, tree, paths)
        result = count(diff) > 0
        finalize(diff)
    catch err
        result = true
    finally
        finalize(tree)
    end
    return result
end

""" git diff --name-only --diff-filter=<filter> <branch1> <branch2> """
function diff_files(repo::GitRepo, branch1::AbstractString, branch2::AbstractString;
                    filter::Set{Cint}=Set([Consts.DELTA_ADDED, Consts.DELTA_MODIFIED, Consts.DELTA_DELETED]))
    b1_id = revparseid(repo, branch1*"^{tree}")
    b2_id = revparseid(repo, branch2*"^{tree}")
    tree1 = get(GitTree, repo, b1_id)
    tree2 = get(GitTree, repo, b2_id)
    files = AbstractString[]
    try
        diff = diff_tree(repo, tree1, tree2)
        for i in 1:count(diff)
            delta = diff[i]
            delta === nothing && break
            if delta.status in filter
                push!(files, unsafe_string(delta.new_file.path))
            end
        end
        finalize(diff)
    finally
        finalize(tree1)
        finalize(tree2)
    end
    return files
end

function is_ancestor_of(a::AbstractString, b::AbstractString, repo::GitRepo)
    A = revparseid(repo, a)
    merge_base(repo, a, b) == A
end

function set_remote_url(repo::GitRepo, url::AbstractString; remote::AbstractString="origin")
    with(GitConfig, repo) do cfg
        set!(cfg, "remote.$remote.url", url)

        m = match(GITHUB_REGEX,url)
        if m !== nothing
            push = "git@github.com:$(m.captures[1]).git"
            if push != url
                set!(cfg, "remote.$remote.pushurl", push)
            end
        end
    end
end

function set_remote_url(path::AbstractString, url::AbstractString; remote::AbstractString="origin")
    with(GitRepo, path) do repo
        set_remote_url(repo, url, remote=remote)
    end
end

function make_payload{P<:AbstractCredentials}(payload::Nullable{P})
    Ref{Nullable{AbstractCredentials}}(payload)
end

""" git fetch [<url>|<repository>] [<refspecs>]"""
function fetch{T<:AbstractString, P<:AbstractCredentials}(repo::GitRepo;
                                  remote::AbstractString="origin",
                                  remoteurl::AbstractString="",
                                  refspecs::Vector{T}=AbstractString[],
                                  payload::Nullable{P}=Nullable{AbstractCredentials}())
    rmt = if isempty(remoteurl)
        get(GitRemote, repo, remote)
    else
        GitRemoteAnon(repo, remoteurl)
    end
    try
        payload = make_payload(payload)
        fo = FetchOptions(callbacks=RemoteCallbacks(credentials_cb(), payload))
        fetch(rmt, refspecs, msg="from $(url(rmt))", options = fo)
    finally
        finalize(rmt)
    end
end

""" git push [<url>|<repository>] [<refspecs>]"""
function push{T<:AbstractString, P<:AbstractCredentials}(repo::GitRepo;
              remote::AbstractString="origin",
              remoteurl::AbstractString="",
              refspecs::Vector{T}=AbstractString[],
              force::Bool=false,
              payload::Nullable{P}=Nullable{AbstractCredentials}())
    rmt = if isempty(remoteurl)
        get(GitRemote, repo, remote)
    else
        GitRemoteAnon(repo, remoteurl)
    end
    try
        payload = make_payload(payload)
        push_opts=PushOptions(callbacks=RemoteCallbacks(credentials_cb(), payload))
        push(rmt, refspecs, force=force, options=push_opts)
    finally
        finalize(rmt)
    end
end

""" git branch """
function branch(repo::GitRepo)
    head_ref = head(repo)
    try
        branch(head_ref)
    finally
        finalize(head_ref)
    end
end

""" git checkout [-b|-B] <branch> [<start-point>] [--track <remote>/<branch>] """
function branch!(repo::GitRepo, branch_name::AbstractString,
                 commit::AbstractString = ""; # start point
                 track::AbstractString  = "", # track remote branch
                 force::Bool=false,           # force branch creation
                 set_head::Bool=true)         # set as head reference on exit
    # try to lookup branch first
    branch_ref = force ? nothing : lookup_branch(repo, branch_name)
    if branch_ref === nothing
        branch_rmt_ref = isempty(track) ? nothing : lookup_branch(repo, "$track/$branch_name", true)
        # if commit is empty get head commit oid
        commit_id = if isempty(commit)
            if branch_rmt_ref === nothing
                with(head(repo)) do head_ref
                    with(peel(GitCommit, head_ref)) do hrc
                        Oid(hrc)
                    end
                end
            else
                tmpcmt = with(peel(GitCommit, branch_rmt_ref)) do hrc
                    Oid(hrc)
                end
                finalize(branch_rmt_ref)
                tmpcmt
            end
        else
            Oid(commit)
        end
        iszero(commit_id) && return
        cmt =  get(GitCommit, repo, commit_id)
        new_branch_ref = nothing
        try
            new_branch_ref = create_branch(repo, branch_name, cmt, force=force)
        finally
            finalize(cmt)
            new_branch_ref === nothing && throw(GitError(Error.Object, Error.ERROR, "cannot create branch `$branch_name` with `$commit_id`"))
            branch_ref = new_branch_ref
        end
    end
    try
        #TODO: what if branch tracks other then "origin" remote
        if !isempty(track) # setup tracking
            try
                with(GitConfig, repo) do cfg
                    set!(cfg, "branch.$branch_name.remote", Consts.REMOTE_ORIGIN)
                    set!(cfg, "branch.$branch_name.merge", name(branch_ref))
                end
            catch
                warn("Please provide remote tracking for branch '$branch_name' in '$(path(repo))'")
            end
        end

        if set_head
            # checkout selected branch
            with(peel(GitTree, branch_ref)) do btree
                checkout_tree(repo, btree)
            end

            # switch head to the branch
            head!(repo, branch_ref)
        end
    finally
        finalize(branch_ref)
    end
    return
end

""" git checkout [-f] --detach <commit> """
function checkout!(repo::GitRepo, commit::AbstractString = "";
                  force::Bool = true)
    # nothing to do
    isempty(commit) && return

    # grab head name
    head_name = Consts.HEAD_FILE
    try
        with(head(repo)) do head_ref
            head_name = shortname(head_ref)
            # if it is HEAD use short OID instead
            if head_name == Consts.HEAD_FILE
                head_name = string(Oid(head_ref))
            end
        end
    end

    # search for commit to get a commit object
    obj = get(GitAnyObject, repo, Oid(commit))
    obj === nothing && return
    try
        peeled = peel(obj, Consts.OBJ_COMMIT)
        peeled === nothing && return
        opts = force ? CheckoutOptions(checkout_strategy = Consts.CHECKOUT_FORCE) :
                       CheckoutOptions()
        try
            # detach commit
            obj_oid = Oid(peeled)
            ref = GitReference(repo, obj_oid, force=force,
                msg="libgit2.checkout: moving from $head_name to $(string(obj_oid))")
            finalize(ref)

            # checkout commit
            checkout_tree(repo, peeled, options = opts)
        finally
            finalize(peeled)
        end
    finally
        finalize(obj)
    end
end

""" git clone [-b <branch>] [--bare] <url> <dir> """
function clone{P<:AbstractCredentials}(repo_url::AbstractString, repo_path::AbstractString;
               branch::AbstractString="",
               isbare::Bool = false,
               remote_cb::Ptr{Void} = C_NULL,
               payload::Nullable{P}=Nullable{AbstractCredentials}())
    # setup clone options
    lbranch = Base.cconvert(Cstring, branch)
    payload = make_payload(payload)
    fetch_opts=FetchOptions(callbacks = RemoteCallbacks(credentials_cb(), payload))
    clone_opts = CloneOptions(
                bare = Cint(isbare),
                checkout_branch = isempty(lbranch) ? Cstring(C_NULL) : Base.unsafe_convert(Cstring, lbranch),
                fetch_opts=fetch_opts,
                remote_cb = remote_cb
            )
    return clone(repo_url, repo_path, clone_opts)
end

""" git reset [<committish>] [--] <pathspecs>... """
function reset!(repo::GitRepo, committish::AbstractString, pathspecs::AbstractString...)
    obj = revparse(repo, !isempty(committish) ? committish : Consts.HEAD_FILE)
    # do not remove entries in the index matching the provided pathspecs with empty target commit tree
    obj === nothing && throw(GitError(Error.Object, Error.ERROR, "`$committish` not found"))
    try
        reset!(repo, Nullable(obj), pathspecs...)
    finally
        finalize(obj)
    end
end

""" git reset [--soft | --mixed | --hard] <commit> """
function reset!(repo::GitRepo, commit::Oid, mode::Cint = Consts.RESET_MIXED)
    obj = get(GitAnyObject, repo, commit)
    # object must exist for reset
    obj === nothing && throw(GitError(Error.Object, Error.ERROR, "Commit `$(string(commit))` object not found"))
    try
        reset!(repo, obj, mode)
    finally
        finalize(obj)
    end
end

""" git cat-file <commit> """
function cat{T<:GitObject}(repo::GitRepo, ::Type{T}, object::AbstractString)
    obj_id = revparseid(repo, object)
    iszero(obj_id) && return nothing

    obj = get(T, repo, obj_id)
    if isa(obj, GitBlob)
        return unsafe_string(convert(Ptr{UInt8}, content(obj)))
    else
        return nothing
    end
end

""" git rev-list --count <commit1> <commit2> """
function revcount(repo::GitRepo, fst::AbstractString, snd::AbstractString)
    fst_id = revparseid(repo, fst)
    snd_id = revparseid(repo, snd)
    base_id = merge_base(repo, string(fst_id), string(snd_id))
    fc = with(GitRevWalker(repo)) do walker
        count((i,r)->i!=base_id, walker, oid=fst_id, by=Consts.SORT_TOPOLOGICAL)
    end
    sc = with(GitRevWalker(repo)) do walker
        count((i,r)->i!=base_id, walker, oid=snd_id, by=Consts.SORT_TOPOLOGICAL)
    end
    return (fc-1, sc-1)
end

""" git merge [--ff-only] [<committish> | FETCH_HEAD] """
function merge!(repo::GitRepo;
                committish::AbstractString = "",
                branch::AbstractString = "",
                fastforward::Bool = false,
                merge_opts::MergeOptions = MergeOptions(),
                checkout_opts::CheckoutOptions = CheckoutOptions())
    # merge into head branch
    with(head(repo)) do head_ref
        upst_anns = if !isempty(committish) # merge committish into HEAD
            if committish == Consts.FETCH_HEAD # merge FETCH_HEAD
                fheads = fetchheads(repo)
                filter!(fh->fh.ismerge, fheads)
                if isempty(fheads)
                    throw(GitError(Error.Merge, Error.ERROR,
                                   "There is no fetch reference for this branch."))
                end
                map(fh->GitAnnotated(repo,fh), fheads)
            else # merge commitish
                [GitAnnotated(repo, committish)]
            end
        else
            if !isempty(branch) # merge provided branch into HEAD
                with(GitReference(repo, branch)) do brn_ref
                    [GitAnnotated(repo, brn_ref)]
                end
            else # try to get tracking remote branch for the head
                if !isattached(repo)
                    throw(GitError(Error.Merge, Error.ERROR,
                                   "Repository HEAD is detached. Remote tracking branch cannot be used."))
                end
                with(upstream(head_ref)) do tr_brn_ref
                    if tr_brn_ref === nothing
                        throw(GitError(Error.Merge, Error.ERROR,
                                       "There is no tracking information for the current branch."))
                    end
                    [GitAnnotated(repo, tr_brn_ref)]
                end
            end
        end

        try
            merge!(repo, upst_anns, fastforward,
                   merge_opts=merge_opts,
                   checkout_opts=checkout_opts)
        finally
            map(finalize, upst_anns)
        end
    end
end

""" git rebase --merge [--onto <newbase>] [<upstream>] """
function rebase!(repo::GitRepo, upstream::AbstractString="", newbase::AbstractString="")
    with(head(repo)) do head_ref
        head_ann = GitAnnotated(repo, head_ref)
        upst_ann  = if isempty(upstream)
            with(upstream(head_ref)) do brn_ref
                if brn_ref === nothing
                    throw(GitError(Error.Rebase, Error.ERROR,
                                   "There is no tracking information for the current branch."))
                end
                GitAnnotated(repo, brn_ref)
            end
        else
            GitAnnotated(repo, upstream)
        end
        try
            sig = default_signature(repo)
            try
                rbs = GitRebase(repo, head_ann, upst_ann)
                try
                    while (rbs_op = next(rbs)) !== nothing
                        commit(rbs, sig)
                    end
                    finish(rbs, sig)
                catch err
                    abort(rbs)
                finally
                    finalize(rbs)
                end
            finally
                #!isnull(onto_ann) && finalize(get(onto_ann))
                finalize(sig)
            end
        finally
            finalize(upst_ann)
            finalize(head_ann)
        end
    end
end

""" Returns all commit authors """
function authors(repo::GitRepo)
    return with(GitRevWalker(repo)) do walker
        map((oid,repo)->with(get(GitCommit, repo, oid)) do cmt
                            author(cmt)::Signature
                        end,
            walker) #, by = Consts.SORT_TIME)
    end
end

function snapshot(repo::GitRepo)
    head = Oid(repo, Consts.HEAD_FILE)
    index = with(GitIndex, repo) do idx; write_tree!(idx) end
    work = try
        with(GitIndex, repo) do idx
            if length(readdir(path(repo))) > 1
                add!(idx, ".")
                write!(idx)
            end
            write_tree!(idx)
        end
    finally
        # restore index
        with(GitIndex, repo) do idx
            read_tree!(idx, index)
            write!(idx)
        end
    end
    State(head, index, work)
end

function restore(s::State, repo::GitRepo)
    reset!(repo, Consts.HEAD_FILE, "*")  # unstage everything
    with(GitIndex, repo) do idx
        read_tree!(idx, s.work)            # move work tree to index
        opts = CheckoutOptions(
                checkout_strategy = Consts.CHECKOUT_FORCE |     # check the index out to work
                                    Consts.CHECKOUT_REMOVE_UNTRACKED) # remove everything else
        checkout_index(repo, Nullable(idx), options = opts)

        read_tree!(idx, s.index)  # restore index
    end
    reset!(repo, s.head, Consts.RESET_SOFT) # restore head
end

function transact(f::Function, repo::GitRepo)
    state = snapshot(repo)
    try f(repo) catch
        restore(state, repo)
        rethrow()
    finally
        finalize(repo)
    end
end

function set_ssl_cert_locations(cert_loc)
    cert_file = isfile(cert_loc) ? cert_loc : Cstring(C_NULL)
    cert_dir  = isdir(cert_loc) ? cert_loc : Cstring(C_NULL)
    cert_file == C_NULL && cert_dir == C_NULL && return
    # TODO FIX https://github.com/libgit2/libgit2/pull/3935#issuecomment-253910017
    #ccall((:git_libgit2_opts, :libgit2), Cint,
    #      (Cint, Cstring, Cstring),
    #      Cint(Consts.SET_SSL_CERT_LOCATIONS), cert_file, cert_dir)
    ENV["SSL_CERT_FILE"] = cert_file
    ENV["SSL_CERT_DIR"] = cert_dir
end

function __init__()
    # Look for OpenSSL env variable for CA bundle (linux only)
    # windows and macOS use the OS native security backends
    old_ssl_cert_dir = Base.get(ENV, "SSL_CERT_DIR", nothing)
    old_ssl_cert_file = Base.get(ENV, "SSL_CERT_FILE", nothing)
    @static if is_linux()
        cert_loc = if "SSL_CERT_DIR" in keys(ENV)
            ENV["SSL_CERT_DIR"]
        elseif "SSL_CERT_FILE" in keys(ENV)
            ENV["SSL_CERT_FILE"]
        else
            # If we have a bundled ca cert file, point libgit2 at that so SSL connections work.
            abspath(ccall(:jl_get_julia_home, Any, ()),Base.DATAROOTDIR,"julia","cert.pem")
        end
        set_ssl_cert_locations(cert_loc)
    end

    err = ccall((:git_libgit2_init, :libgit2), Cint, ())
    err > 0 || throw(ErrorException("error initializing LibGit2 module"))
    atexit() do
        ccall((:git_libgit2_shutdown, :libgit2), Cint, ())
    end

    @static if is_linux()
        if old_ssl_cert_dir != Base.get(ENV, "SSL_CERT_DIR", "")
            if old_ssl_cert_dir === nothing
                delete!(ENV, "SSL_CERT_DIR")
            else
                ENV["SSL_CERT_DIR"] = old_ssl_cert_dir
            end
        end
        if old_ssl_cert_file != Base.get(ENV, "SSL_CERT_FILE", "")
            if old_ssl_cert_file === nothing
                delete!(ENV, "SSL_CERT_FILE")
            else
                ENV["SSL_CERT_FILE"] = old_ssl_cert_file
            end
        end
    end
end

end # module

max_stars_repo_name: habemus-papadum/julia

max_issues_repo_name: habemus-papadum/julia

max_forks_repo_name: habemus-papadum/julia


text:
using Test
using ParameterEstimocean
using ParameterEstimocean.Utils: prettyvector

@test prettyvector([0, 1, 2]) == "[0, 1, 2]"
@test prettyvector(collect(0:20)) == "[0, 1, 2 … 20, 19, 18] (21 elements)"


text:
using ImageFiltering

@testset "watershed" begin

    img = zeros(100, 100)
    img[25:75, 25:75] .= 1
    img[26:74, 26:74] .= 0

    markers = zeros(Int, size(img))
    markers[1, 1] = 1
    markers[50, 50] = 2

    @testset "classic watershed" begin
        # with classic watershed, we expect the top left label should spread around
        # the center label.
        result = watershed(img, markers)

        @test length(result.segment_labels) == 2
        @test result.segment_labels == collect(1:2)
        @test all(label->(label==result.image_indexmap[50,50]), result.image_indexmap[26:74,26:74])
    end

    @testset "masked watershed" begin
        mask = trues(size(img))
        mask[60:70, 60:70] .= false

        result = watershed(img, markers, mask=mask)
        labels = labels_map(result)

        # where the mask is false, no label should be assigned
        @test sum(labels[.~ mask]) == 0

        result = watershed(img, markers, compactness=10.0, mask=mask)
        labels = labels_map(result)

        # where the mask is false, no label should be assigned
        @test sum(labels[.~ mask]) == 0
    end

    @testset "compact watershed" begin
        result = watershed(img, markers, compactness=10.0)
        labels = labels_map(result)

        # since this is using the compact algorithm with a high value for
        # compactness, the boundary between labels 1 and 2 should occur halfway
        # between the two markers
        @test sum(labels .== 1) == sum(1:50)
    end

    @testset "h-minima transform" begin
        img = ones(15, 15)
        #minima of depth 0.2
        img[3:5, 3:5] .= 0.9
        img[4,4] = 0.8
        #minima of depth 0.7
        img[9:13, 9:13] .= 0.8
        img[10:12, 10:12] .= 0.7
        img[11,11] = 0.3

        out = hmin_transform(img, 0.25)

        @test findlocalminima(img) == [CartesianIndex(4, 4), CartesianIndex(11, 11)]
        @test findlocalminima(out) == [CartesianIndex(11, 11)]
    end
end


text:
module SamplingMethods

# Types
abstract type SamplingMethod <: Function end
abstract type SamplingMethodState <: Function end

const CuSVTuple = Tuple{
    AnyCuArray, AnyCuArray, AnyCuArray, AnyCuArray, AnyCuArray, NamedTuple, Real
}

const Snapshot{T} = NamedTuple{
    (:positions, :velocities, :masses, :forces, :ids, :box, :dt), T
}

include("abf.jl")

end  # module SamplingMethods


text:
# using Base.Test
# using MIToS.MSA

print("""

Test shuffle_...
================
""")

let aln = read(joinpath(pwd(), "data", "PF09645_full.fasta.gz"), FASTA), copy_aln_msa = copy(aln.msa),
  gaps = (aln .== GAP), lcol = mean(aln .== Residue('L'), 1), lseq = mean(aln .== Residue('L'), 2)

  shuffle_residues_columnwise!(aln)
  @test (aln .== GAP) == gaps
  @test aln.msa != copy_aln_msa
  @test lcol == mean(aln .== Residue('L'), 1)
  @test lseq != mean(aln .== Residue('L'), 2)
end

let aln = read(joinpath(pwd(), "data", "PF09645_full.fasta.gz"), FASTA), copy_aln_msa = copy(aln.msa),
  gaps = (aln .== GAP), lcol = mean(aln .== Residue('L'), 1), lseq = mean(aln .== Residue('L'), 2)

  shuffle_residues_sequencewise!(aln)
  @test (aln .== GAP) == gaps
  @test aln.msa != copy_aln_msa
  @test lcol != mean(aln .== Residue('L'), 1)
  @test lseq == mean(aln .== Residue('L'), 2)
end

let aln = read(joinpath(pwd(), "data", "PF09645_full.fasta.gz"), FASTA), copy_aln_msa = copy(aln.msa),
  gaps = (aln .== GAP), lcol = mean(aln .== Residue('L'), 1), lseq = mean(aln .== Residue('L'), 2)

  shuffle_columnwise!(aln)
  @test (aln .== GAP) != gaps
  @test aln.msa != copy_aln_msa
  @test lcol == mean(aln .== Residue('L'), 1)
  @test lseq != mean(aln .== Residue('L'), 2)
end

let aln = read(joinpath(pwd(), "data", "PF09645_full.fasta.gz"), FASTA), copy_aln_msa = copy(aln.msa),
  gaps = (aln .== GAP), lcol = mean(aln .== Residue('L'), 1), lseq = mean(aln .== Residue('L'), 2)

  shuffle_sequencewise!(aln)
  @test (aln .== GAP) != gaps
  @test aln.msa != copy_aln_msa
  @test lcol != mean(aln .== Residue('L'), 1)
  @test lseq == mean(aln .== Residue('L'), 2)
end


text:
julia> meandegrees([350, 10])
0.0

julia> meandegrees([90, 180, 270, 360]])
0.0

julia> meandegrees([10, 20, 30]])
19.999999999999996


text:
using PolarFact
using Test, LinearAlgebra
using Random

Random.seed!(1103)

include("test_newton.jl")
include("test_halley.jl")
include("test_svd.jl")
include("test_hybrid.jl")
include("test_f32.jl") # test Float32 case


text:
#!/usr/local/bin/julia
# last tested Julia version: 1.6.1
# Morris-Lecar model on a 2D lattice
# FvW 03/2018

using NPZ
using PyCall
using PyPlot
using Statistics
using VideoIO
@pyimport matplotlib.animation as anim

function ml2d(N, T, t0, dt, s, D, gL, VL, gCa, VCa, gK, VK, C, I0,
			  V0, V1, V2, V3, V4, phi, stim, blocks)
    C1 = 1.0/C
    # initialize Morris-Lecar system
    v = V0*ones(Float64,N,N)
    w = 0.015*ones(Float64,N,N)
    dv = zeros(Float64,N,N)
    dw = zeros(Float64,N,N)
    s_sqrt_dt = s*sqrt(dt)
    X = zeros(Float64,T,N,N)
    # stimulation protocol
    I = zeros(Float64,t0+T,N,N)
    for st in stim
        t_on, t_off = st[1]
        x0, x1 = st[2]
        y0, y1 = st[3]
        I[t0+t_on:t0+t_off, x0:x1, y0:y1] .= I0
    end

    # iterate
    for t in range(1, stop=t0+T, step=1)
        (t%100 == 0) && print("    t = ", t, "/", t0+T, "\r")
        # ML equations
        dv = C1.*(-gL.*(v.-VL) -gCa.*0.5.*(1.0 .+ tanh.((v.-V1)/V2)).*(v.-VCa)
                 -gK.*w.*(v.-VK) .+ I[t,:,:]) .+ D.*L(v)
        dw = phi*cosh.((v.-V3)/(2*V4)) .* (0.5*(1.0 .+ tanh.((v.-V3)/V4))-w)
        # Ito stochastic integration
        v += (dv*dt + s_sqrt_dt*randn(N,N))
        w += (dw*dt)
        # dead block(s):
        for bl in blocks
            v[bl[1][1]:bl[1][2], bl[2][1]:bl[2][2]] .= 0.0
            w[bl[1][1]:bl[1][2], bl[2][1]:bl[2][2]] .= 0.0
        end
        (t > t0) && (X[t-t0,:,:] = v)
    end
    println("\n")
    return X
end

function animate_pyplot(fname, data)
    """
    Animate 3D array as .mp4 using PyPlot, save as `fname`
    array dimensions:
        1: time
        2, 3: space
    """
    nt, nx, ny = size(data)
    vmin = minimum(data)
    vmax = maximum(data)
    # setup animation image
    fig = figure(figsize=(6,6))
    axis("off")
    t = imshow(data[1,:,:], origin="lower", cmap=ColorMap("gray"),
			   vmin=vmin, vmax=vmax)
    tight_layout()
    # frame generator
    println("[+] animate")
    function animate(i)
        (i%100 == 0) && print("    t = ", i, "/", nt, "\r")
        t.set_data(data[i+1,:,:])
    end
    # create animation
    ani = anim.FuncAnimation(fig, animate, frames=nt, interval=10)
    println("\n")
    # save animation
    ani[:save](fname, bitrate=-1,
               extra_args=["-vcodec", "libx264", "-pix_fmt", "yuv420p"])
    show()
end

function animate_video(fname, data)
    """
    Animate 3D array as .mp4 using VideoIO, save as `fname`
    array dimensions:
        1: time
        2, 3: space
    """
    # BW
    y = UInt8.(round.(255*(data .- minimum(data))/(maximum(data)-minimum(data))))
    # BW inverted
    #y = UInt8.(round.(255 .- 255*(data .- minimum(data))/(maximum(data)-minimum(data))))
    encoder_options = (color_range=2, crf=0, preset="medium")
    framerate=30
    T = size(data,1)
    open_video_out(fname, y[1,end:-1:1,:], framerate=framerate,
                   encoder_options=encoder_options) do writer
        for i in range(2,stop=T,step=1)
            write(writer, y[i,end:-1:1,:])
        end
    end
end

function L(x)
    # Laplace operator
    # periodic boundary conditions
    xU = circshift(x, [-1 0])
    xD = circshift(x, [1 0])
    xL = circshift(x, [0 -1])
    xR = circshift(x, [0 1])
    Lx = xU + xD + xL + xR - 4x
    # non-periodic boundary conditions
    Lx[1,:] .= 0.0
    Lx[end,:] .= 0.0
    Lx[:,1] .= 0.0
    Lx[:,end] .= 0.0
    return Lx
end

function main()
    println("Morris-Lecar (ML) lattice model\n")
    N = 128
    T = 15000
    t0 = 1000
	dt = 0.01
	s = 2.0
	D = 0.8
	# Morris-Lecar parameters
    gL = 2.0
    VL = -60.0
    gCa = 5.5 # 4.0
    VCa = 120.0
    gK = 8.0
    VK = -84.0
    C = 5.0 # 20.0
    I = 120.0
    V0 = -60.0
    V1 = -1.2
    V2 = 18.0
    V3 = 12.0 # 2.0
    V4 = 17.4 # 30.0
    phi = 1.0/15.0 # 0.04
    println("[+] Lattice size N: ", N)
    println("[+] Time steps T: ", T)
    println("[+] Warm-up steps t0: ", t0)
	println("[+] Integration time step dt: ", dt)
	println("[+] Noise intensity s: ", s)
	println("[+] Diffusion coefficient D: ", D)
    println("[+] ML parameter gL: ", gL)
    println("[+] ML parameter VL: ", VL)
    println("[+] ML parameter gCa: ", gCa)
    println("[+] ML parameter VCa: ", VCa)
    println("[+] ML parameter gK: ", gK)
    println("[+] ML parameter VK: ", VK)
    println("[+] ML parameter C: ", C)
    println("[+] ML parameter V0: ", V0)
    println("[+] ML parameter V1: ", V1)
    println("[+] ML parameter V2: ", V2)
    println("[+] ML parameter V3: ", V3)
    println("[+] ML parameter V4: ", V4)
    println("[+] ML parameter phi: ", phi)
	println("[+] Stimulation current I: ", I)

    # stim protocol, array of elements [[t0,t1], [x0,x1], [y0,y1]]
    # double spiral 1
    #stim = [ [[50,250], [1,5], [1,5]],
    #         [[5800,6200], [20,25], [1,35]] ]
    # double spiral 2
    #stim = [ [[50,250], [1,5], [1,5]],
    #         [[5800,6200], [18,23], [5,26]] ]
    # spiral wave 1
    #stim = [ [[50,250], [1,5], [1,5]],
    #         [[6300,6800], [30,35], [5,20]] ]
    # spiral wave 2
    stim = [ [[50,250], [1,5], [1,5]],
             [[6400,6900], [30,35], [5,20]] ]
    #stim = []

    # dead blocks, array of elementy [[x0,x1], [y0,y1]]
    blocks = [ [[1,15], [5,10]] ]
    #blocks = []

    # run simulation
    data = ml2d(N, T, t0, dt, s, D, gL, VL, gCa, VCa, gK, VK, C, I,
				V0, V1, V2, V3, V4, phi, stim, blocks)
    println("[+] Data dimensions: ", size(data))

    # plot mean voltage
    m = mean(reshape(data, (T,N*N)), dims=2)
    plot(m, "-k"); show()

    # save data
    I_str = rpad(I, 4, '0') # stim. current amplitude as 4-char string
    s_str = rpad(s, 4, '0') # noise as 4-char string
    D_str = rpad(D, 4, '0') # diffusion coefficient as 4-char string
    #fname1 = string("ml2d_I_", I_str, "_s_", s_str, "_D_", D_str, ".npy")
    #npzwrite(fname1, data)
    #println("[+] Data saved as: ", fname1)

    # video
    fname2 = string("ml2d_I_", I_str, "_s_", s_str, "_D_", D_str, ".mp4")
    #animate_pyplot(fname2, data) # slow
    animate_video(fname2, data) # fast
    println("[+] Data saved as: ", fname2)
end

main()


text:
# -*- coding: utf-8 -*-
# ---
# jupyter:
#   jupytext:
#     formats: ipynb,jl:light
#     text_representation:
#       extension: .jl
#       format_name: light
#       format_version: '1.4'
#       jupytext_version: 1.2.4
#   kernelspec:
#     display_name: Julia 1.3.1
#     language: julia
#     name: julia-1.3
# ---

# + {"slideshow": {"slide_type": "slide"}, "cell_type": "markdown"}
# # Ocean transports and maps
#
# Transports in the climate system are often represented as gridded vector fields (e.g. on a `C-grid`) and integrated across `grid edge paths` that e.g. (1) connect location pairs or (2) tracks latitude circles. For more detail, please refer to [Forget et al, 2015](https://doi.org/10.5194/gmd-8-3071-2015) (incl. appendices)

# + {"slideshow": {"slide_type": "slide"}, "cell_type": "markdown"}
# ### Read grid & transports from file
#
# 1. pre-requisites
# 2. read variables

# +
#]add MITgcmTools#master

# + {"slideshow": {"slide_type": "subslide"}}
using MeshArrays, Plots, Statistics, MITgcmTools
include(joinpath(dirname(pathof(MeshArrays)),"../examples/Plots.jl"))
include("helper_functions.jl");

get_grid_if_needed()
Γ=GridLoad(GridSpec("LatLonCap","../inputs/GRID_LLC90/"))
(Tx,Ty,τx,τy,η)=trsp_read("LatLonCap","../inputs/GRID_LLC90/");

# + {"slideshow": {"slide_type": "slide"}, "cell_type": "markdown"}
# ### Integrate transport across latitude lines
#
# 1. `LatitudeCircles` computes `grid edge path`s that track latitude circles
# 2. `ThroughFlow` integrates transports accross the specified `grid edge path`s
# 3. Plot integrated meridional transport in `Sverdrup` units

# + {"slideshow": {"slide_type": "subslide"}}
uv=Dict("U"=>Tx,"V"=>Ty,"dimensions"=>["x","y"])
L=-89.0:89.0; LC=LatitudeCircles(L,Γ)

T=Array{Float64,1}(undef,length(LC))
[T[i]=1e-6*ThroughFlow(uv,LC[i],Γ) for i=1:length(LC)]

plot(L,T,xlabel="latitude",ylabel="Northward Transport (in Sv=10^6 m³/s)",label="ECCO")

# + {"slideshow": {"slide_type": "slide"}, "cell_type": "markdown"}
# ## Transport Directions
#
# 1. `u,v` are oriented in the `Eastward,Northward` directions
# 2. `uC,vC` are oriented along the `x,y` directions of each subdomain

# + {"slideshow": {"slide_type": "-"}}
u,v,uC,vC=rotate_uv(uv,Γ);

# + {"slideshow": {"slide_type": "subslide"}, "cell_style": "split"}
heatmap(u,clims=(-20.0,20.0),title="East-ward")
#heatmap(v,clims=(-20.0,20.0),title="North-ward")

# + {"slideshow": {"slide_type": "fragment"}, "cell_style": "split"}
heatmap(uC,clims=(-20.0,20.0),title="x-ward")
#heatmap(vC,clims=(-20.0,20.0),title="y-ward")

# + {"slideshow": {"slide_type": "slide"}, "cell_type": "markdown"}
# ### Global Maps
#
# 1. interpolate `u,v` to a `1/2 x 1/2` degree grid for plotting
# 2. map out the Eastward,Northward transport components

# + {"slideshow": {"slide_type": "-"}}
uI,vI,lon,lat=interp_uv(u,v);
# + {"slideshow": {"slide_type": "subslide"}, "cell_style": "split"}
heatmap(lon,lat,uI,clims=(-20.0,20.0),
    title="Eastward transport (in Sv / cell)")
# + {"slideshow": {"slide_type": "fragment"}, "cell_style": "split"}
heatmap(lon,lat,vI,clims=(-20.0,20.0),
    title="Northward transport (in Sv / cell)")

max_stars_repo_name: abelolivelanna/GlobalOceanNotebooks

max_issues_repo_name: abelolivelanna/GlobalOceanNotebooks

max_forks_repo_name: abelolivelanna/GlobalOceanNotebooks


text:
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#
# Description
# ==============================================================================
#
#   Tests related to the beta angle analysis.
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #

# File: ./src/beta_angle.jl
# ==============================================================================

# Function: beta_angle
# ------------------------------------------------------------------------------

@testset "Function beta_angle" begin
  jd₀ = SatelliteAnalysis.date_to_jd(2021, 1, 1, 0, 0, 0)

  β = beta_angle(
      jd₀,
      7130.982e3,
      0.001111,
      deg2rad(98.405),
      SatelliteAnalysis.compute_RAAN_lt(jd₀, 22.5),
      5
  )

  β_expected = [
      0.43291445352510594
      0.43481148866086317
      0.43668726746800246
      0.4385397145910499
      0.44036676126144947
  ]

  @test β ≈ β_expected

  β = beta_angle(
      jd₀,
      7130.982e3,
      0.001111,
      deg2rad(98.405),
      SatelliteAnalysis.compute_RAAN_lt(jd₀, 22.5) + π,
      5
  )

  β_expected = [
      -0.31068187881262466
      -0.31295402147649054
      -0.31524372027870284
      -0.3175489312409039
      -0.3198676044177273
  ]

  @test β ≈ β_expected
end

@testset "Function beta_angle (errors)" begin
  jd₀ = SatelliteAnalysis.date_to_jd(2021, 1, 1, 0, 0, 0)

  @test_throws Exception beta_angle(
      jd₀,
      7130.982e3,
      0.001111,
      deg2rad(98.405),
      SatelliteAnalysis.compute_RAAN_lt(jd₀, 22.5),
      0
  )

  @test_throws Exception beta_angle(
      jd₀,
      7130.982e3,
      0.001111,
      deg2rad(98.405),
      SatelliteAnalysis.compute_RAAN_lt(jd₀, 22.5),
      -5
  )
end


text:
@testset "swing" begin
    tree = ParseNewick("((I,J)A,(K,(M,N)L)B,(C,((O,P)D,E)F)G)H;")

    l = tree_length(tree)

    swing!(tree)
    t2 = swing(tree)
    l1 = tree_length(tree)
    l2 = tree_length(t2)
    # take care of numerical instabilities
    @test l == l1
    @test l == l2
end

@testset "slide" begin
    tree = ParseNewick("((I,J)A,(K,(M,N)L)B,(C,((O,P)D,E)F)G)H;")

    l = tree_length(tree)

    slide!(tree)
    l2 = tree_length(tree)

    t2 = slide(tree)
    l3 = tree_length(t2)
    # take care of numerical instabilities
    @test l == round(l2)
    @test l == round(l3)
end

@testset "NNI" begin
    tree = ParseNewick("((I,J)A,(K,(M,N)L)B,(C,(O,P)E)G)H;")

    tr2 = deepcopy(tree)

    res = NNI!(tr2, 10, true)

    l = tree_length(tree)
    l2 = tree_length(tr2)

    @test l == l2
    @test MCPhyloTree.RF_int(tree, tr2)[1] == 2
end

@testset "move" begin
    tree = ParseNewick("((B:5,A:5)C:9,(D:5,E:5)F:5)G:5;")
    B = find_by_name(tree,"B")
    A = find_by_name(tree,"A")
    MCPhyloTree.move!(B,A,0.9)
    @test B.inc_length == 9.0
    @test A.inc_length == 1.0
end

@testset "change_edge_length" begin
    tree = ParseNewick("((B:1,A:1)C:1,(D:1,E:1)F:1)G:1;")
    change_edge_length!(tree)
    po_list = [x.inc_length for x in post_order(tree)]
    @test sum(po_list) != 7.0
end

@testset "reroot" begin
    tree = ParseNewick("((B:1,A:1)C:1,(D:1,E:1)F:1)G:1;")
    newtree = MCPhyloTree.reroot(tree,"C")
    @test newick(newtree) == "(B:1.0,A:1.0,((D:1.0,E:1.0)F:1.0)G:1.0)C:1.0;"
end

@testset "SPR" begin
    tree2 = ParseNewick("(A,B,(C,(D,E)F)G)H;")
    subtree = find_by_name(tree2,"F")
    target = find_by_name(tree2,"A")
    nutree = MCPhyloTree.perform_spr(tree2,subtree,target)
    @test newick(nutree) == "(B:1.0,C:1.0,((D:1.0,E:1.0)F:1.0,A:1.0)G:1.0)H:1.0;"
end