text:
using Clang
using Clang.Generators
using Clang.LibClang.Clang_jll
using Test
using Clang.Generators: strip_comment_markers

include("rewriter.jl")

@testset "Generators" begin
    INCLUDE_DIR = normpath(Clang_jll.artifact_dir, "include")
    CLANG_C_DIR = joinpath(INCLUDE_DIR, "clang-c")

    options = load_options(joinpath(@__DIR__, "test.toml"))

    # add compiler flags
    args = get_default_args()
    push!(args, "-I$INCLUDE_DIR")

    # search top-level headers
    headers = detect_headers(CLANG_C_DIR, args)

    # add extra definition
    @add_def time_t AbstractJuliaSIT JuliaCtime_t Ctime_t

    # create context
    ctx = create_context(headers, args, options)

    # build without printing so we can do custom rewriting
    build!(ctx, BUILDSTAGE_NO_PRINTING)

    rewrite!(ctx.dag)

    # print
    @test_logs (:info, "Done!") match_mode=:any build!(ctx, BUILDSTAGE_PRINTING_ONLY)
end

@testset "Comments" begin
    @test strip_comment_markers("/* abc */") == ["abc "]
    @test strip_comment_markers("/** abc */") == ["abc "]
    @test strip_comment_markers("/*< abc */") == ["abc "]
    @test strip_comment_markers("/// hello") == ["hello"]
    @test strip_comment_markers("/**\n * line1\n * line2\n */") == ["line1", "line2"]
    @test strip_comment_markers("/*!\n * line1\n * line2\n */") == ["line1", "line2"]
    @test strip_comment_markers("    /// line1\n    /// line2") == ["line1", "line2"]
    @test strip_comment_markers("//! line1\n//! line2") == ["line1", "line2"]
    @test strip_comment_markers("//! line1") == ["line1"]
    @test strip_comment_markers("//< line1") == ["line1"]
end

@testset "Resolve dependency" begin
    args = get_default_args()
    headers = joinpath(@__DIR__, "include", "dependency.h")
    options = Dict("general" => Dict{String,Any}(
            "output_file_path" => joinpath(@__DIR__, "LibDependency.jl")))
    ctx = create_context(headers, args, options)
    build!(ctx)
    @test include("LibDependency.jl") isa Any
end

@testset "Issue 320" begin
    args = get_default_args()
    dir = joinpath(@__DIR__, "sys")
    push!(args, "-isystem$dir")
    headers = [joinpath(@__DIR__, "include", "test.h")]
    ctx = create_context(headers, args)
    @add_def stat
    @test build!(ctx, BUILDSTAGE_NO_PRINTING) isa Any
end

@testset "Escape anonymous name with var\"\"" begin
    args = get_default_args()
    headers = joinpath(@__DIR__, "include", "escape-with-var.h")
    options = Dict("general" => Dict{String,Any}(
            "output_file_path" => joinpath(@__DIR__, "LibEscapeWithVar.jl")))
    ctx = create_context(headers, args, options)
    build!(ctx)
    @test include("LibEscapeWithVar.jl") isa Any
end

@testset "Issue 307" begin
    args = get_default_args()
    dir = joinpath(@__DIR__, "sys")
    push!(args, "-isystem$dir")
    headers = joinpath(@__DIR__, "include", "struct-in-union.h")
    ctx = create_context(headers, args)
    @test build!(ctx, BUILDSTAGE_NO_PRINTING) isa Any

    headers = joinpath(@__DIR__, "include", "nested-struct.h")
    ctx = create_context(headers, args)
    @test build!(ctx, BUILDSTAGE_NO_PRINTING) isa Any

    headers = joinpath(@__DIR__, "include", "nested-declaration.h")
    ctx = create_context(headers, args)
    @test_broken try
        build!(ctx, BUILDSTAGE_NO_PRINTING) isa Any
        true
    catch
        false
    end
end

@testset "Issue 327" begin
    tu = parse_header(Index(), joinpath(@__DIR__, "include/void-type.h"))
    root = Clang.getTranslationUnitCursor(tu)
    func = children(root)[]
    ret_type = Clang.getCursorResultType(func)
    @test ret_type isa CLVoid
end

@testset "Issue 355" begin
    tu = parse_header(Index(), joinpath(@__DIR__, "include/return-funcptr.h"))
    root = Clang.getTranslationUnitCursor(tu)
    func = children(root)[3]
    @test length(get_function_args(func)) == 1
end

@testset "macros" begin
    ctx = create_context(joinpath(@__DIR__, "include/macro.h"), get_default_args())
    @test_logs (:info, "Done!") match_mode = :any build!(ctx)
end

@testset "#368" begin
    ctx = create_context(joinpath(@__DIR__, "include/union-in-struct.h"),
                         get_default_args())
    @test_logs (:info, "Done!") match_mode = :any build!(ctx)
    @test ctx.dag.nodes[end].id == :A
    @test ctx.dag.nodes[end].type isa Generators.StructLayout
end

@testset "Issue 376" begin
    ctx = create_context(joinpath(@__DIR__, "include/macro-dependency.h"), get_default_args())
    @test build!(ctx) isa Any
end


text:
# Test tracing functionalities

module TracingTest
    using ContextTracking
    using Test

    @ctx function foo()
        @memo x = 1
        me = "not visible downstream"
        bar()
        @test context().path == [:foo]
    end

    @ctx function bar()
        @memo y = 2
        z = 3
        @memo z
        baz()
        @test context().path == [:foo, :bar]
    end

    @ctx function baz()
        c = context()
        @test c.path == [:foo, :bar, :baz]
        @test c.data[:x] == 1
        @test c.data[:y] == 2
        @test c.data[:z] == 3
        @test !haskey(c.data, :me)
    end

    __init__() = foo()
end

ext: jl


text:
using Test, LMPTools
using Dates
using GeographicLib
using LongwaveModePropagator
const LMP = LongwaveModePropagator

const GROUND_DATA = LMPTools.GROUND_DATA
const GROUND_LAT = GROUND_DATA["lat"]
const GROUND_LON = GROUND_DATA["lon"]

function bfieldlatlons(tx, rx)
    dists = 0:100e3:range(tx, rx)
    line = GeodesicLine(tx, rx)

    lats = Vector{Float64}(undef, length(dists))
    lons = similar(lats)
    for i in eachindex(dists)
        for (i, d) in enumerate(dists)
            lo, la, _ = forward(line, d)
            lats[i] = la
            lons[i] = lo
        end
    end
    return dists, lats, lons
end

@testset "LMPTools" begin
    tx = TRANSMITTER[:NAA]
    rx = Receiver("Boulder", 40.01, -105.244, 0.0, VerticalDipole())

    # Extended functions
    gl1 = GeodesicLine(tx, rx)
    gl2 = GeodesicLine(tx.longitude, tx.latitude; lon2=rx.longitude, lat2=rx.latitude)
    for f in fieldnames(GeodesicLine)
        # testing with ≈ because some values are ~1e-315
        @test getfield(gl1, f) ≈ getfield(gl2, f)
    end

    # https://www.cqsrg.org/tools/GCDistance/
    @test range(tx, rx) ≈ 3142.0e3 atol=100

    # Ground
    @test extrema(GROUND_DATA["lat"]) == (-90, 90)
    @test extrema(GROUND_DATA["lon"]) == (-180, 180)

    @test LMPTools.searchsortednearest(GROUND_DATA["lat"], 78) ==
        findfirst(isequal(78), GROUND_DATA["lat"])
    @test LMPTools.searchsortednearest(GROUND_DATA["lat"], 78.1) ==
        findfirst(isequal(78), GROUND_DATA["lat"])
    @test LMPTools.searchsortednearest(GROUND_DATA["lat"], 78.4) ==
        findfirst(isequal(78.5), GROUND_DATA["lat"])
    @test LMPTools.searchsortednearest(GROUND_DATA["lat"], -78) ==
        findfirst(isequal(-78), GROUND_DATA["lat"])
    @test LMPTools.searchsortednearest(GROUND_DATA["lat"], -78.1) ==
        findfirst(isequal(-78), GROUND_DATA["lat"])
    @test LMPTools.searchsortednearest(GROUND_DATA["lat"], -78.4) ==
        findfirst(isequal(-78.5), GROUND_DATA["lat"])

    for la in GROUND_LAT, lo in GROUND_LON
        @test get_ground("epsilonmap", la, lo) == get_epsilon(la, lo)
        @test get_ground("sigmamap", la, lo) == get_sigma(la, lo)
        @test get_ground("codemap", la, lo) == get_groundcode(la, lo)

        @test get_epsilon(la, lo) ≈ GROUND[get_groundcode(la, lo)].ϵᵣ
        @test get_sigma(la, lo) ≈ GROUND[get_groundcode(la, lo)].σ atol=1e-8
    end

    grounds, distances = groundsegments(tx, rx)
    @test minimum(diff(distances)) >= 20e3
    @test length(grounds) != length(unique(grounds))
    @test length(grounds) == length(distances)

    grounds, distances = groundsegments(tx, rx; require_permittivity_change=true)
    @test length(grounds) == length(unique(grounds))
    @test length(grounds) == length(distances)

    # Magnetic field

    az, _, dist, _ = inverse(tx.longitude, tx.latitude, rx.longitude, rx.latitude)

    # Based on https://www.ngdc.noaa.gov/geomag/calculators/magcalc.shtml#igrfwmm
    declination = -16.0378
    calcb = BField(50.6701e-6, deg2rad(67.7737), deg2rad(az-declination))

    bfield = igrf(az, tx.latitude, tx.longitude, 2020)
    @test all(isapprox(getfield(bfield,f), getfield(calcb,f); rtol=1e-4) for f in fieldnames(BField))

    bfield2 = only(igrf(tx, rx, 2020, 0.0))
    @test all(isapprox(getfield(bfield2,f), getfield(bfield,f)) for f in fieldnames(BField))

    # chaos does not match igrf, but they should be reasonably close
    bfield3 = chaos(az, tx.latitude, tx.longitude, 2020)
    @test all(isapprox(getfield(bfield3,f), getfield(bfield,f); rtol=0.1) for f in fieldnames(BField))

    ranges, lats, lons = bfieldlatlons(tx, rx)
    chaos_bfields = chaos(tx, rx, 2020, ranges)
    igrf_bfields = igrf(tx, rx, 2020, ranges)

    @test chaos_bfields == chaos(az, lats, lons, 2020)
    for i in 1:length(ranges)
        @test chaos_bfields[i] == chaos(az, lats[i], lons[i], 2020)
        @test all(isapprox(getfield(chaos_bfields[i],f), getfield(igrf_bfields[i],f); rtol=1) for f in fieldnames(BField))
    end

    # Try loading a different CHAOS mat file
    load_CHAOS_matfile(joinpath(LMPTools.project_path("data"), "CHAOS-7.8.mat"))
    bfield4 = chaos(az, tx.latitude, tx.longitude, 2020)
    @test all(isapprox(getfield(bfield4,f), getfield(bfield,f); rtol=0.1) for f in fieldnames(BField))

    # Zenith angle

    dt = DateTime(2021, 2, 1, 12, 45)
    lon, lat = forward(gl1, 1000e3)
    @test zenithangle(lat, lon, dt) == zenithangle(lat, lon, 2021, 2, 1, 12.75)

    # Check out https://ssd.jpl.nasa.gov/horizons.cgi for comparison numbers
    @test zenithangle(lat, lon, dt) ≈ 90 - 0.3106 atol=0.01
    @test zenithangle(lat, lon, 2021, 2, 1, 16) ≈ 90 - 24.9418 atol=0.01
    @test zenithangle(lat, lon, 2021, 2, 1, 21.25) ≈ 90 - 10.5537 atol=0.01

    lon, lat = -110, 62
    @test zenithangle(lat, lon, 2021, 2, 1, 9.5) ≈ 90 - -40.5603 atol=0.01
    @test zenithangle(lat, lon, 2021, 2, 1, 20) ≈ 90 - 10.9651 atol=0.01
    @test zenithangle(lat, lon, 2021, 9, 15, 5) ≈ 90 - -20.1376 atol=0.01
    @test zenithangle(lat, lon, 2021, 9, 15, 19.25) ≈ 90 - 30.7122 atol=0.01

    @test isday(zenithangle(lat, lon, 2021, 2, 1, 16)) == true  # day
    @test isday(zenithangle(lat, lon, 2021, 2, 1, 22)) == true
    @test isday(zenithangle(lat, lon, 2021, 2, 2, 1)) == false  # night
    @test isday(zenithangle(lat, lon, 2021, 2, 2, 6)) == false
    @test isday(zenithangle(lat, lon, 2021, 2, 2, 8)) == false
    @test isday(zenithangle(lat, lon, 2021, 2, 2, 11)) == false

    # Ionospheres

    dt = DateTime(2020, 3, 1, 2)
    lats = 15:89
    lons = -160:-60
    szas = [zenithangle(la, lo, dt) for la in lats, lo in lons]

    # Compare number and array-specialized form
    x = [flatlinearterminator(s) for s in szas]
    hprimes, betas = getindex.(x,1), getindex.(x,2)
    hprimes_arr, betas_arr = flatlinearterminator(szas)
    @test hprimes == hprimes_arr
    @test betas == betas_arr

    h, b = smoothterminator(50)
    @test h ≈ 74  # hp_day
    @test b ≈ 0.3  # b_day

    h, b = smoothterminator(130)
    @test h ≈ 86  # hp_night
    @test b ≈ 0.5  # b_night

    h, b = smoothterminator(95)
    @test 74 < h < 86
    @test 0.3 < b < 0.5

    # Sanity tests on fourierperturbation
    coeff = zeros(5)
    pert = fourierperturbation.(szas, (coeff,))
    @test iszero(pert)

    coeff = [1, 0, 0, 0, 0]
    pert = fourierperturbation.(szas, (coeff,))
    @test all(isequal(1), pert)
end


text:
@testset "Static graphs" begin
    g = @inferred(CompleteDiGraph(5))
    @test nv(g) == 5 && ne(g) == 20
    g = @inferred(CompleteGraph(5))
    @test nv(g) == 5 && ne(g) == 10

    g = @inferred(CompleteBipartiteGraph(5, 8))
    @test nv(g) == 13 && ne(g) == 40

    g = @inferred(StarDiGraph(5))
    @test nv(g) == 5 && ne(g) == 4
    g = @inferred(StarGraph(5))
    @test nv(g) == 5 && ne(g) == 4
    g = @inferred(StarGraph(1))
    @test nv(g) == 1 && ne(g) == 0

    g = @inferred(PathDiGraph(5))
    @test nv(g) == 5 && ne(g) == 4
    g = @inferred(PathGraph(5))
    @test nv(g) == 5 && ne(g) == 4

    g = @inferred(CycleDiGraph(5))
    @test nv(g) == 5 && ne(g) == 5
    g = @inferred(CycleGraph(5))
    @test nv(g) == 5 && ne(g) == 5

    g = @inferred(WheelDiGraph(5))
    @test nv(g) == 5 && ne(g) == 8
    g = @inferred(WheelGraph(5))
    @test nv(g) == 5 && ne(g) == 8

    g = @inferred(Grid([3, 3, 4]))
    @test nv(g) == 3 * 3 * 4
    @test ne(g) == 75
    @test Δ(g) == 6
    @test δ(g) == 3

    g = @inferred(Grid([3, 3, 4], periodic=true))
    @test nv(g) == 3 * 3 * 4
    @test ne(g) == 108
    @test Δ(g) == 6
    @test δ(g) == 6

    g = @inferred(CliqueGraph(3, 5))
    @test nv(g) == 15 && ne(g) == 20
    @test g[1:3] == CompleteGraph(3)

    g = @inferred(crosspath(3, BinaryTree(2)))
    # f = Vector{Vector{Int}}[[2 3 4];
    # [1 5];
    # [1 6];
    # [1 5 6 7];
    # [2 4 8];
    # [3 4 9];
    # [4 8 9];
    # [5 7];
    # [6 7]
    # ]
    I = [1, 1, 1, 2, 2, 3, 3, 4, 4, 4, 4, 5, 5, 5, 6, 6, 6, 7, 7, 7, 8, 8, 9, 9]
    J = [2, 3, 4, 1, 5, 1, 6, 1, 5, 6, 7, 2, 4, 8, 3, 4, 9, 4, 8, 9, 5, 7, 6, 7]
    V = ones(Int, length(I))
    Adj = sparse(I, J, V)
    @test Adj == sparse(g)

    g = @inferred(DoubleBinaryTree(3))
    # [[3, 2, 8]
    # [4, 1, 5]
    # [1, 6, 7]
    # [2]
    # [2]
    # [3]
    # [3]
    # [10, 9, 1]
    # [11, 8, 12]
    # [8, 13, 14]
    # [9]
    # [9]
    # [10]
    # [10]]
    I = [1, 1, 1, 2, 2, 2, 3, 3, 3, 4, 5, 6, 7, 8, 8, 8, 9, 9, 9, 10, 10, 10, 11, 12, 13, 14]
    J = [3, 2, 8, 4, 1, 5, 1, 6, 7, 2, 2, 3, 3, 10, 9, 1, 11, 8, 12, 8, 13, 14, 9, 9, 10, 10]
    V = ones(Int, length(I))
    Adj = sparse(I, J, V)
    @test Adj == sparse(g)

    rg3 = @inferred(RoachGraph(3))
    # [3]
    # [4]
    # [1, 5]
    # [2, 6]
    # [3, 7]
    # [4, 8]
    # [9, 8, 5]
    # [10, 7, 6]
    # [11, 10, 7]
    # [8, 9, 12]
    # [9, 12]
    # [10, 11]
    I = [1, 2, 3, 3, 4, 4, 5, 5, 6, 6, 7, 7, 7, 8, 8, 8, 9, 9, 9, 10, 10, 10, 11, 11, 12, 12]
    J = [3, 4, 1, 5, 2, 6, 3, 7, 4, 8, 9, 8, 5, 10, 7, 6, 11, 10, 7, 8, 9, 12, 9, 12, 10, 11]
    V = ones(Int, length(I))
    Adj = sparse(I, J, V)
    @test Adj == sparse(rg3)
end


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

abstract type AbstractTime end

"""
    Period
    Year
    Month
    Week
    Day
    Hour
    Minute
    Second
    Millisecond
    Microsecond
    Nanosecond

`Period` types represent discrete, human representations of time.
"""
abstract type Period     <: AbstractTime end
abstract type DatePeriod <: Period end
abstract type TimePeriod <: Period end

for T in (:Year, :Month, :Week, :Day)
    @eval struct $T <: DatePeriod
        value::Int64
        $T(v::Number) = new(v)
    end
end
for T in (:Hour, :Minute, :Second, :Millisecond, :Microsecond, :Nanosecond)
    @eval struct $T <: TimePeriod
        value::Int64
        $T(v::Number) = new(v)
    end
end

"""
    Year(v)
    Month(v)
    Week(v)
    Day(v)
    Hour(v)
    Minute(v)
    Second(v)
    Millisecond(v)
    Microsecond(v)
    Nanosecond(v)

Construct a `Period` type with the given `v` value. Input must be losslessly convertible
to an [`Int64`](@ref).
"""
Period(v)

"""
    Instant

`Instant` types represent integer-based, machine representations of time as continuous
timelines starting from an epoch.
"""
abstract type Instant <: AbstractTime end

"""
    UTInstant{T}

The `UTInstant` represents a machine timeline based on UT time (1 day = one revolution of
the earth). The `T` is a `Period` parameter that indicates the resolution or precision of
the instant.
"""
struct UTInstant{P<:Period} <: Instant
    periods::P
end

# Convenience default constructors
UTM(x) = UTInstant(Millisecond(x))
UTD(x) = UTInstant(Day(x))

# Calendar types provide rules for interpretating instant
# timelines in human-readable form.
abstract type Calendar <: AbstractTime end

# ISOCalendar implements the ISO 8601 standard (en.wikipedia.org/wiki/ISO_8601)
# Notably based on the proleptic Gregorian calendar
# ISOCalendar provides interpretation rules for UTInstants to civil date and time parts
struct ISOCalendar <: Calendar end

abstract type TimeZone end
struct UTC <: TimeZone end

"""
    TimeType

`TimeType` types wrap `Instant` machine instances to provide human representations of the
machine instant. `Time`, `DateTime` and `Date` are subtypes of `TimeType`.
"""
abstract type TimeType <: AbstractTime end

abstract type AbstractDateTime <: TimeType end

"""
    DateTime

`DateTime` wraps a `UTInstant{Millisecond}` and interprets it according to the proleptic
Gregorian calendar.
"""
struct DateTime <: AbstractDateTime
    instant::UTInstant{Millisecond}
    DateTime(instant::UTInstant{Millisecond}) = new(instant)
end

"""
    Date

`Date` wraps a `UTInstant{Day}` and interprets it according to the proleptic Gregorian calendar.
"""
struct Date <: TimeType
    instant::UTInstant{Day}
    Date(instant::UTInstant{Day}) = new(instant)
end

"""
    Time

`Time` wraps a `Nanosecond` and represents a specific moment in a 24-hour day.
"""
struct Time <: TimeType
    instant::Nanosecond
    Time(instant::Nanosecond) = new(instant)
end

# Convert y,m,d to # of Rata Die days
# Works by shifting the beginning of the year to March 1,
# so a leap day is the very last day of the year
const SHIFTEDMONTHDAYS = (306, 337, 0, 31, 61, 92, 122, 153, 184, 214, 245, 275)
function totaldays(y, m, d)
    # If we're in Jan/Feb, shift the given year back one
    z = m < 3 ? y - 1 : y
    mdays = SHIFTEDMONTHDAYS[m]
    # days + month_days + year_days
    return d + mdays + 365z + fld(z, 4) - fld(z, 100) + fld(z, 400) - 306
end

# If the year is divisible by 4, except for every 100 years, except for every 400 years
isleapyear(y) = ((y % 4 == 0) && (y % 100 != 0)) || (y % 400 == 0)

# Number of days in month
const DAYSINMONTH = (31, 28, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31)
daysinmonth(y,m) = DAYSINMONTH[m] + (m == 2 && isleapyear(y))

### UTILITIES ###

# These are necessary because the type constructors for TimeType subtypes can
# throw, and we want to be able to use tryparse without requiring a try/catch.
# This is made easier by providing a helper function that checks arguments, so
# we can validate arguments in tryparse.

"""
    validargs(::Type{<:TimeType}, args...) -> Union{ArgumentError, Nothing}

Determine whether the given arguments consitute valid inputs for the given type.
Returns either an `ArgumentError`, or [`nothing`](@ref) in case of success.
"""
function validargs end

"""
    argerror([msg]) -> Union{ArgumentError, Nothing}

Return an `ArgumentError` object with the given message,
or [`nothing`](@ref) if no message is provided. For use by `validargs`.
"""
argerror(msg::String) = ArgumentError(msg)
argerror() = nothing

### CONSTRUCTORS ###
# Core constructors
"""
    DateTime(y, [m, d, h, mi, s, ms]) -> DateTime

Construct a `DateTime` type by parts. Arguments must be convertible to [`Int64`](@ref).
"""
function DateTime(y::Int64, m::Int64=1, d::Int64=1,
                  h::Int64=0, mi::Int64=0, s::Int64=0, ms::Int64=0)
    err = validargs(DateTime, y, m, d, h, mi, s, ms)
    err === nothing || throw(err)
    rata = ms + 1000 * (s + 60mi + 3600h + 86400 * totaldays(y, m, d))
    return DateTime(UTM(rata))
end

function validargs(::Type{DateTime}, y::Int64, m::Int64, d::Int64,
                   h::Int64, mi::Int64, s::Int64, ms::Int64)
    0 < m < 13 || return argerror("Month: $m out of range (1:12)")
    0 < d < daysinmonth(y, m) + 1 || return argerror("Day: $d out of range (1:$(daysinmonth(y, m)))")
    -1 < h < 24 || (h == 24 && mi==s==ms==0) ||
        return argerror("Hour: $h out of range (0:23)")
    -1 < mi < 60 || return argerror("Minute: $mi out of range (0:59)")
    -1 < s < 60 || return argerror("Second: $s out of range (0:59)")
    -1 < ms < 1000 || return argerror("Millisecond: $ms out of range (0:999)")
    return argerror()
end

"""
    Date(y, [m, d]) -> Date

Construct a `Date` type by parts. Arguments must be convertible to [`Int64`](@ref).
"""
function Date(y::Int64, m::Int64=1, d::Int64=1)
    err = validargs(Date, y, m, d)
    err === nothing || throw(err)
    return Date(UTD(totaldays(y, m, d)))
end

function validargs(::Type{Date}, y::Int64, m::Int64, d::Int64)
    0 < m < 13 || return argerror("Month: $m out of range (1:12)")
    0 < d < daysinmonth(y, m) + 1 || return argerror("Day: $d out of range (1:$(daysinmonth(y, m)))")
    return argerror()
end

"""
    Time(h, [mi, s, ms, us, ns]) -> Time

Construct a `Time` type by parts. Arguments must be convertible to [`Int64`](@ref).
"""
function Time(h::Int64, mi::Int64=0, s::Int64=0, ms::Int64=0, us::Int64=0, ns::Int64=0)
    err = validargs(Time, h, mi, s, ms, us, ns)
    err === nothing || throw(err)
    return Time(Nanosecond(ns + 1000us + 1000000ms + 1000000000s + 60000000000mi + 3600000000000h))
end

function validargs(::Type{Time}, h::Int64, mi::Int64, s::Int64, ms::Int64, us::Int64, ns::Int64)
    -1 < h < 24 || return argerror("Hour: $h out of range (0:23)")
    -1 < mi < 60 || return argerror("Minute: $mi out of range (0:59)")
    -1 < s < 60 || return argerror("Second: $s out of range (0:59)")
    -1 < ms < 1000 || return argerror("Millisecond: $ms out of range (0:999)")
    -1 < us < 1000 || return argerror("Microsecond: $us out of range (0:999)")
    -1 < ns < 1000 || return argerror("Nanosecond: $ns out of range (0:999)")
    return argerror()
end

# Convenience constructors from Periods
function DateTime(y::Year, m::Month=Month(1), d::Day=Day(1),
                  h::Hour=Hour(0), mi::Minute=Minute(0),
                  s::Second=Second(0), ms::Millisecond=Millisecond(0))
    return DateTime(value(y), value(m), value(d),
                        value(h), value(mi), value(s), value(ms))
end

Date(y::Year, m::Month=Month(1), d::Day=Day(1)) = Date(value(y), value(m), value(d))

function Time(h::Hour, mi::Minute=Minute(0), s::Second=Second(0),
              ms::Millisecond=Millisecond(0),
              us::Microsecond=Microsecond(0), ns::Nanosecond=Nanosecond(0))
    return Time(value(h), value(mi), value(s), value(ms), value(us), value(ns))
end

# To allow any order/combination of Periods

"""
    DateTime(periods::Period...) -> DateTime

Construct a `DateTime` type by `Period` type parts. Arguments may be in any order. DateTime
parts not provided will default to the value of `Dates.default(period)`.
"""
function DateTime(period::Period, periods::Period...)
    y = Year(1); m = Month(1); d = Day(1)
    h = Hour(0); mi = Minute(0); s = Second(0); ms = Millisecond(0)
    for p in (period, periods...)
        isa(p, Year) && (y = p::Year)
        isa(p, Month) && (m = p::Month)
        isa(p, Day) && (d = p::Day)
        isa(p, Hour) && (h = p::Hour)
        isa(p, Minute) && (mi = p::Minute)
        isa(p, Second) && (s = p::Second)
        isa(p, Millisecond) && (ms = p::Millisecond)
    end
    return DateTime(y, m, d, h, mi, s, ms)
end

"""
    Date(period::Period...) -> Date

Construct a `Date` type by `Period` type parts. Arguments may be in any order. `Date` parts
not provided will default to the value of `Dates.default(period)`.
"""
function Date(period::Period, periods::Period...)
    y = Year(1); m = Month(1); d = Day(1)
    for p in (period, periods...)
        isa(p, Year) && (y = p::Year)
        isa(p, Month) && (m = p::Month)
        isa(p, Day) && (d = p::Day)
    end
    return Date(y, m, d)
end

"""
    Time(period::TimePeriod...) -> Time

Construct a `Time` type by `Period` type parts. Arguments may be in any order. `Time` parts
not provided will default to the value of `Dates.default(period)`.
"""
function Time(period::TimePeriod, periods::TimePeriod...)
    h = Hour(0); mi = Minute(0); s = Second(0)
    ms = Millisecond(0); us = Microsecond(0); ns = Nanosecond(0)
    for p in (period, periods...)
        isa(p, Hour) && (h = p::Hour)
        isa(p, Minute) && (mi = p::Minute)
        isa(p, Second) && (s = p::Second)
        isa(p, Millisecond) && (ms = p::Millisecond)
        isa(p, Microsecond) && (us = p::Microsecond)
        isa(p, Nanosecond) && (ns = p::Nanosecond)
    end
    return Time(h, mi, s, ms, us, ns)
end

# Fallback constructors
DateTime(y, m=1, d=1, h=0, mi=0, s=0, ms=0) = DateTime(Int64(y), Int64(m), Int64(d), Int64(h), Int64(mi), Int64(s), Int64(ms))
Date(y, m=1, d=1) = Date(Int64(y), Int64(m), Int64(d))
Time(h, mi=0, s=0, ms=0, us=0, ns=0) = Time(Int64(h), Int64(mi), Int64(s), Int64(ms), Int64(us), Int64(ns))

# Traits, Equality
Base.isfinite(::Union{Type{T}, T}) where {T<:TimeType} = true
calendar(dt::DateTime) = ISOCalendar
calendar(dt::Date) = ISOCalendar

"""
    eps(::DateTime) -> Millisecond
    eps(::Date) -> Day
    eps(::Time) -> Nanosecond

Returns `Millisecond(1)` for `DateTime` values, `Day(1)` for `Date` values, and `Nanosecond(1)` for `Time` values.
"""
Base.eps

Base.eps(dt::DateTime) = Millisecond(1)
Base.eps(dt::Date) = Day(1)
Base.eps(t::Time) = Nanosecond(1)

Base.typemax(::Union{DateTime, Type{DateTime}}) = DateTime(146138512, 12, 31, 23, 59, 59)
Base.typemin(::Union{DateTime, Type{DateTime}}) = DateTime(-146138511, 1, 1, 0, 0, 0)
Base.typemax(::Union{Date, Type{Date}}) = Date(252522163911149, 12, 31)
Base.typemin(::Union{Date, Type{Date}}) = Date(-252522163911150, 1, 1)
Base.typemax(::Union{Time, Type{Time}}) = Time(23, 59, 59, 999, 999, 999)
Base.typemin(::Union{Time, Type{Time}}) = Time(0)
# Date-DateTime promotion, isless, ==
Base.promote_rule(::Type{Date}, x::Type{DateTime}) = DateTime
Base.isless(x::T, y::T) where {T<:TimeType} = isless(value(x), value(y))
Base.isless(x::TimeType, y::TimeType) = isless(promote(x, y)...)
(==)(x::T, y::T) where {T<:TimeType} = (==)(value(x), value(y))
function ==(a::Time, b::Time)
    return hour(a) == hour(b) && minute(a) == minute(b) &&
        second(a) == second(b) && millisecond(a) == millisecond(b) &&
        microsecond(a) == microsecond(b) && nanosecond(a) == nanosecond(b)
end
(==)(x::TimeType, y::TimeType) = (===)(promote(x, y)...)
Base.hash(x::Time, h::UInt) =
    hash(hour(x), hash(minute(x), hash(second(x),
        hash(millisecond(x), hash(microsecond(x), hash(nanosecond(x), h))))))

import Base: sleep, Timer, timedwait
sleep(time::Period) = sleep(toms(time) / 1000)
Timer(time::Period; interval::Period = Second(0)) =
    Timer(toms(time) / 1000, interval = toms(interval) / 1000)
timedwait(testcb::Function, time::Period) = timedwait(testcb, toms(time) / 1000)

Base.OrderStyle(::Type{<:AbstractTime}) = Base.Ordered()
Base.ArithmeticStyle(::Type{<:AbstractTime}) = Base.ArithmeticWraps()

max_stars_repo_name: skumagai/julia

max_issues_repo_name: skumagai/julia

max_forks_repo_name: skumagai/julia


text:
"""

Taken from Vazza et. al. (2016): doi:10.1093/mnras/stw584

They assume Caprioli+Spitkovsky (2014) efficiency is roughly half of Kang&Ryu (2013)

We account for angle dependent efficiency independently.
"""
const vazza_factor = 1//2

"""
        CS14(X_cr::T=0.05, η_max::T=0.5*0.2055) where T

Efficiency model by Caprioli&Spitkovsky 2014 efficiency.

## Values
- `X_cr`: P_cr / P_th defined in model for re-acceleration. Basis for interpolation between acceleration and re-acceleration efficiency.
- `η_max`: Maximum efficiency defined in the model
"""
struct CS14{T} <: AbstractShockAccelerationEfficiency
    X_cr::T
    η_max::T

    CS14(X_cr::T = 0.05, η_max::T = 0.5 * 0.2055) where {T} = new{T}(X_cr, η_max)
end

"""
    η_Ms_acc(η_model::CS14, M::T) where T

Efficiency from Caprioli&Spitkovsky 2015, doi: 10.1088/0004-637x/783/2/91
Same simplified approach as Vazza+16 -> is roughly half the efficiency of Kang&Ryu 2013.
"""
function η_Ms_acc(η_model::CS14, M::Real)
    return vazza_factor * η_Ms_acc(KR13(), M)
end

"""
    η_Ms_reacc(η_model::CS14, M::T) where T

Efficiency from Caprioli&Spitkovsky 2015, doi: 10.1088/0004-637x/783/2/91
Same simplified approach as Vazza+16 -> is roughly half the efficiency of Kang&Ryu 2013.
"""
function η_Ms_reacc(η_model::CS14, M::Real)
    return vazza_factor * η_Ms_reacc(KR13(), M)
end


text:
using Expronicon
using Documenter

DocMeta.setdocmeta!(Expronicon, :DocTestSetup, :(using Expronicon); recursive=true)

makedocs(;
    modules=[Expronicon],
    authors="Roger-luo <rogerluo.rl18@gmail.com> and contributors",
    repo="https://github.com/Roger-luo/Expronicon.jl/blob/{commit}{path}#{line}",
    sitename="Expronicon.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://Roger-luo.github.io/Expronicon.jl",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
        "Types" => "types.md",
        "Analysis" => "analysis.md",
        "Transform" => "transform.md",
        "CodeGen" => "codegen.md",
        "Printings" => "printings.md",
    ],
)

deploydocs(;
    repo="github.com/Roger-luo/Expronicon.jl",
)


text:
using ProteinEnsembles
using Test

testfile(path::AbstractString...) = joinpath(dirname(@__FILE__), "test_files", path...)

# Whether to run the parameterisation tests
# This is to satisfy the online auto-builds
# Set to true locally to run the tests
run_param_test = false

if run_param_test
    # Optional argument is path to TMscore executable
    if isempty(ARGS)
        test_tmscore_path = "TMscore"
    else
        test_tmscore_path = ARGS[1]
    end
    @assert tmscorepathvalid(test_tmscore_path) "Not a valid TMscore path: \"$test_tmscore_path\""
    println("TMscore executable path is taken as \"$test_tmscore_path\"")
end

# Temp directory which is removed at the end
temp_dir = "$(tempdir())/$(defaults["out_dir"])"

include("test_types.jl")
include("test_atoms.jl")
include("test_io.jl")
include("test_interactions.jl")
include("test_generate.jl")
include("test_align.jl")
include("test_pca.jl")
include("test_perturb.jl")
include("test_pipeline.jl")


text:
@recipe function f(grid::Grid{T, N, S}) where {T, N, S}
    #only plot grid point densities for axes with more than one tick
    iaxs = findall(ax -> length(ax) > 1, grid.axes)
    layout --> length(iaxs)
    seriestype --> :stephist
    legend --> false

    for (i,iax) in enumerate(iaxs)
        @series begin
            subplot := i
            bins --> div(length(grid[iax]), 2)
            xguide --> "Grid point density - Axis $(iax)"
            grid[iax]
        end
    end
end


text:
using LowRankModels
import StatsBase: sample, WeightVec

# test quadratic loss

## generate data
srand(1);
m,n,k = 300,300,3;
kfit = k+1
# variance of measurement
sigmasq = .1

# coordinates of covariates
X_real = randn(m,k)
# directions of observations
Y_real = randn(k,n)

XY = X_real*Y_real;
A = zeros(Int, (m, n))
logistic(x) = 1/(1+exp(-x))
for i=1:m
	for j=1:n
		A[i,j] = logistic(XY[i,j]) >= rand() ? 1 : -1
	end
end

# and the model
losses = LogisticLoss()
rx, ry = QuadReg(.1), QuadReg(.1);
glrm = GLRM(A,losses,rx,ry,kfit)
#scale=false, offset=false, X=randn(kfit,m), Y=randn(kfit,n));

# fit w/o initialization
@time X,Y,ch = fit!(glrm);
XYh = X'*Y;
println("After fitting, parameters differ from true parameters by $(vecnorm(XY - XYh)/sqrt(prod(size(XY)))) in RMSE")
A_imputed = impute(glrm)
println("After fitting, $(sum(A_imputed .!= A) / prod(size(A))*100)\% of imputed entries are wrong")
println("After fitting, imputed entries are off by $(sum(abs(A_imputed - A)) / prod(size(A))*100)\% on average")
println("(Picking randomly, 50\% of entries would be wrong.)\n")

# initialize
init_svd!(glrm)
XYh = glrm.X' * glrm.Y
println("After initialization with the svd, parameters differ from true parameters by $(vecnorm(XY - XYh)/sqrt(prod(size(XY)))) in RMSE")
A_imputed = impute(glrm)
println("After fitting, $(sum(A_imputed .!= A) / prod(size(A))*100)\% of imputed entries are wrong")
println("After fitting, imputed entries are off by $(sum(abs(A_imputed - A)) / prod(size(A))*100)\% on average")
println("(Picking randomly, 50\% of entries would be wrong.)\n")

# fit w/ initialization
@time X,Y,ch = fit!(glrm);
XYh = X'*Y;
println("After fitting, parameters differ from true parameters by $(vecnorm(XY - XYh)/sqrt(prod(size(XY)))) in RMSE")
A_imputed = impute(glrm)
println("After fitting, $(sum(A_imputed .!= A) / prod(size(A))*100)\% of imputed entries are wrong")
println("After fitting, imputed entries are off by $(sum(abs(A_imputed - A)) / prod(size(A))*100)\% on average")
println("(Picking randomly, 50\% of entries would be wrong.)\n")


text:
# Function for plotting the performance of a specific kernel in order to illustrate
# performance difference between PMM and DRAM
using nGraph
using nGraph.Flux
using nGraph.NNlib

# This is the implementation of the convolution kernel for profiling.
# We do some addition tricks to make sure intermediate tensors exist.
function __kernel(x, w)
    x += ones(eltype(x), size(x))
    w += ones(eltype(x), size(w))
    z = NNlib.conv(x, w; pad = (1, 1, 1, 1))
    return log(z)
end

function profile_kernel()
    # Input data
    x = rand(Float32, 48, 48, 128, 128)
    # Input weights
    w = rand(Float32, 3, 3, 128, 128)

    # Create an Actualizer and do the profiling.
    f = () -> AutoTM.Actualizer(__kernel, x, w)
    kernel_profile(f)
end

# Plot the kernel performance
function plot_kernel_performance(;
        file = "plot.pdf",
        width = "15cm",
        height = "5cm",
    )

    # Step 1: Get the kernel from the cache
    cache = AutoTM.Profiler.CPUKernelCache(AutoTM.Experiments.CPU_CACHE)
    configs = (AutoTM.Utils.DRAM, AutoTM.Utils.PMEM)

    expected_input_sizes = (
        (48, 48, 128, 128),
        (3, 3, 128, 128),
    )

    times = Float64[]
    for i in Iterators.product(configs, configs)
        for o in configs
            config = AutoTM.Utils.IOConfig(wrap(i), wrap(o))
            for (k,v) in cache.cache
                params = k[1]
                io_config = k[2]
                if config == io_config &&
                        params.description == "Convolution" &&
                        params.input_sizes == expected_input_sizes
                    push!(times, v)
                end
            end
        end
    end

    # Normalize times to a multple of the "all dram" case.
    times .= times ./ first(times)

    # Now that we've found the kernel times, we generate the plot!!
    xticks = [
        "DRAM\\\\DRAM\\\\DRAM",
        "DRAM\\\\DRAM\\\\\\textbf{PMM}",
        "DRAM\\\\\\textbf{PMM}\\\\DRAM",
        "DRAM\\\\\\textbf{PMM}\\\\\\textbf{PMM}",

        "\\textbf{PMM}\\\\DRAM\\\\DRAM",
        "\\textbf{PMM}\\\\DRAM\\\\\\textbf{PMM}",
        "\\textbf{PMM}\\\\\\textbf{PMM}\\\\DRAM",
        "\\textbf{PMM}\\\\\\textbf{PMM}\\\\\\textbf{PMM}",
    ]

    plot = Plot(Coordinates((i, t) for (i,t) in enumerate(times)))

    plt = TikzDocument()
        push!(plt, """
    \\pgfplotsset{
        width=$width,
        height=$height
    }
    """)

    tikz = TikzPicture()
        axs = @pgf Axis(
        {
            ybar,
            enlarge_x_limits=0.10,
            nodes_near_coords_align={vertical},
            ylabel="Performance relative\\\\to all IO in DRAM",
            ymajorgrids,
            ymin=0,
            ylabel_style={
                align = "center",
            },

            xtick="data",
            xticklabels = xticks,
            xticklabel_style={
                align = "center",
            },
            yticklabel_style={
                "/pgf/number format/fixed",
                "/pgf/number format/precision=5",
            },
            bar_width="20pt",
            clip=false,
        },
        plot,
        """
        \\node[align=right] at (axis cs: -0.1, -0.78) {
        \\textbf{Data In:}\\\\
        \\textbf{Weight:}\\\\
        \\textbf{Data Out:}\\\\
        };
        """,
    )

    push!(tikz, axs)
    push!(plt, tikz)
    pgfsave(file, plt)
    return plt
end

max_stars_repo_name: darchr/AutoTM

max_issues_repo_name: darchr/AutoTM

max_forks_repo_name: darchr/AutoTM


text:
const K_SP = ' '
const K_exclam = '!'
const K_quotedbl = '\"'
const K_numbersign = '#'
const K_dollar = '\$'
const K_percent = '%'
const K_ampersand = '&'
const K_apostrophe = '\''
const K_parentleft = '('
const K_parentright = ')'
const K_asterisk = '*'
const K_plus = '+'
const K_comma = ','
const K_minus = '-'
const K_period = '.'
const K_slash = '/'
const K_0 = '0'
const K_1 = '1'
const K_2 = '2'
const K_3 = '3'
const K_4 = '4'
const K_5 = '5'
const K_6 = '6'
const K_7 = '7'
const K_8 = '8'
const K_9 = '9'
const K_colon = ':'
const K_semicolon = ';'
const K_less = '<'
const K_equal = '='
const K_greater = '>'
const K_question = '?'
const K_at = '@'
const K_A = 'A'
const K_B = 'B'
const K_C = 'C'
const K_D = 'D'
const K_E = 'E'
const K_F = 'F'
const K_G = 'G'
const K_H = 'H'
const K_I = 'I'
const K_J = 'J'
const K_K = 'K'
const K_L = 'L'
const K_M = 'M'
const K_N = 'N'
const K_O = 'O'
const K_P = 'P'
const K_Q = 'Q'
const K_R = 'R'
const K_S = 'S'
const K_T = 'T'
const K_U = 'U'
const K_V = 'V'
const K_W = 'W'
const K_X = 'X'
const K_Y = 'Y'
const K_Z = 'Z'
const K_bracketleft = '['
const K_backslash = '\\'
const K_bracketright = ']'
const K_circum = '^'
const K_underscore = '_'
const K_grave = '`'
const K_a = 'a'
const K_b = 'b'
const K_c = 'c'
const K_d = 'd'
const K_e = 'e'
const K_f = 'f'
const K_g = 'g'
const K_h = 'h'
const K_i = 'i'
const K_j = 'j'
const K_k = 'k'
const K_l = 'l'
const K_m = 'm'
const K_n = 'n'
const K_o = 'o'
const K_p = 'p'
const K_q = 'q'
const K_r = 'r'
const K_s = 's'
const K_t = 't'
const K_u = 'u'
const K_v = 'v'
const K_w = 'w'
const K_x = 'x'
const K_y = 'y'
const K_z = 'z'
const K_braceleft = '{'
const K_bar = '|'
const K_braceright = '}'
const K_tilde = '~'
# Skipping MacroDefinition: iup_isprint(_c)((_c)>31&&(_c)<127)
const K_BS = '\b'
const K_TAB = '\t'
const K_LF = '\n'
const K_CR = '\r'
const K_quoteleft = K_grave
const K_quoteright = K_apostrophe
#const isxkey = iup_isXkey
# Skipping MacroDefinition: iup_isXkey(_c)((_c)>=128)
const K_PAUSE = 0xFF13
const K_ESC = 0xFF1B
const K_HOME = 0xFF50
const K_LEFT = 0xFF51
const K_UP = 0xFF52
const K_RIGHT = 0xFF53
const K_DOWN = 0xFF54
const K_PGUP = 0xFF55
const K_PGDN = 0xFF56
const K_END = 0xFF57
const K_MIDDLE = 0xFF0B
const K_Print = 0xFF61
const K_INS = 0xFF63
const K_Menu = 0xFF67
const K_DEL = 0xFFFF
const K_F1 = 0xFFBE
const K_F2 = 0xFFBF
const K_F3 = 0xFFC0
const K_F4 = 0xFFC1
const K_F5 = 0xFFC2
const K_F6 = 0xFFC3
const K_F7 = 0xFFC4
const K_F8 = 0xFFC5
const K_F9 = 0xFFC6
const K_F10 = 0xFFC7
const K_F11 = 0xFFC8
const K_F12 = 0xFFC9
const K_LSHIFT = 0xFFE1
const K_RSHIFT = 0xFFE2
const K_LCTRL = 0xFFE3
const K_RCTRL = 0xFFE4
const K_LALT = 0xFFE9
const K_RALT = 0xFFEA
const K_NUM = 0xFF7F
const K_SCROLL = 0xFF14
const K_CAPS = 0xFFE5
const K_ccedilla = 0x00E7
const K_Ccedilla = 0x00C7
const K_acute = 0x00B4
const K_diaeresis = 0x00A8
# Skipping MacroDefinition: iup_isShiftXkey(_c)((_c)&0x10000000)
# Skipping MacroDefinition: iup_isCtrlXkey(_c)((_c)&0x20000000)
# Skipping MacroDefinition: iup_isAltXkey(_c)((_c)&0x40000000)
# Skipping MacroDefinition: iup_isSysXkey(_c)((_c)&0x80000000)
# Skipping MacroDefinition: iup_XkeyBase(_c)((_c)&0x0FFFFFFF)
# Skipping MacroDefinition: iup_XkeyShift(_c)((_c)|0x10000000)
# Skipping MacroDefinition: iup_XkeyCtrl(_c)((_c)|0x20000000)
# Skipping MacroDefinition: iup_XkeyAlt(_c)((_c)|0x40000000)
# Skipping MacroDefinition: iup_XkeySys(_c)((_c)|0x80000000)

# TO BE CONTINUED 			JL
macro iup_XkeyShift(c)
	return $c | 0x10000000
end

#=
const K_sHOME = iup_XkeyShift
const K_sUP = iup_XkeyShift
const K_sPGUP = iup_XkeyShift
const K_sLEFT = iup_XkeyShift
const K_sMIDDLE = iup_XkeyShift
const K_sRIGHT = iup_XkeyShift
const K_sEND = iup_XkeyShift
const K_sDOWN = iup_XkeyShift
const K_sPGDN = iup_XkeyShift
const K_sINS = iup_XkeyShift
const K_sDEL = iup_XkeyShift
const K_sSP = iup_XkeyShift
const K_sTAB = iup_XkeyShift
const K_sCR = iup_XkeyShift
const K_sBS = iup_XkeyShift
const K_sPAUSE = iup_XkeyShift
const K_sESC = iup_XkeyShift
const K_sF1 = iup_XkeyShift
const K_sF2 = iup_XkeyShift
const K_sF3 = iup_XkeyShift
const K_sF4 = iup_XkeyShift
const K_sF5 = iup_XkeyShift
const K_sF6 = iup_XkeyShift
const K_sF7 = iup_XkeyShift
const K_sF8 = iup_XkeyShift
const K_sF9 = iup_XkeyShift
const K_sF10 = iup_XkeyShift
const K_sF11 = iup_XkeyShift
const K_sF12 = iup_XkeyShift
const K_sPrint = iup_XkeyShift
const K_sMenu = iup_XkeyShift
const K_cHOME = iup_XkeyCtrl
const K_cUP = iup_XkeyCtrl
const K_cPGUP = iup_XkeyCtrl
const K_cLEFT = iup_XkeyCtrl
const K_cMIDDLE = iup_XkeyCtrl
const K_cRIGHT = iup_XkeyCtrl
const K_cEND = iup_XkeyCtrl
const K_cDOWN = iup_XkeyCtrl
const K_cPGDN = iup_XkeyCtrl
const K_cINS = iup_XkeyCtrl
const K_cDEL = iup_XkeyCtrl
const K_cSP = iup_XkeyCtrl
const K_cTAB = iup_XkeyCtrl
const K_cCR = iup_XkeyCtrl
const K_cBS = iup_XkeyCtrl
const K_cPAUSE = iup_XkeyCtrl
const K_cESC = iup_XkeyCtrl
const K_cCcedilla = iup_XkeyCtrl
const K_cF1 = iup_XkeyCtrl
const K_cF2 = iup_XkeyCtrl
const K_cF3 = iup_XkeyCtrl
const K_cF4 = iup_XkeyCtrl
const K_cF5 = iup_XkeyCtrl
const K_cF6 = iup_XkeyCtrl
const K_cF7 = iup_XkeyCtrl
const K_cF8 = iup_XkeyCtrl
const K_cF9 = iup_XkeyCtrl
const K_cF10 = iup_XkeyCtrl
const K_cF11 = iup_XkeyCtrl
const K_cF12 = iup_XkeyCtrl
const K_cPrint = iup_XkeyCtrl
const K_cMenu = iup_XkeyCtrl
const K_mHOME = iup_XkeyAlt
const K_mUP = iup_XkeyAlt
const K_mPGUP = iup_XkeyAlt
const K_mLEFT = iup_XkeyAlt
const K_mMIDDLE = iup_XkeyAlt
const K_mRIGHT = iup_XkeyAlt
const K_mEND = iup_XkeyAlt
const K_mDOWN = iup_XkeyAlt
const K_mPGDN = iup_XkeyAlt
const K_mINS = iup_XkeyAlt
const K_mDEL = iup_XkeyAlt
const K_mSP = iup_XkeyAlt
const K_mTAB = iup_XkeyAlt
const K_mCR = iup_XkeyAlt
const K_mBS = iup_XkeyAlt
const K_mPAUSE = iup_XkeyAlt
const K_mESC = iup_XkeyAlt
const K_mCcedilla = iup_XkeyAlt
const K_mF1 = iup_XkeyAlt
const K_mF2 = iup_XkeyAlt
const K_mF3 = iup_XkeyAlt
const K_mF4 = iup_XkeyAlt
const K_mF5 = iup_XkeyAlt
const K_mF6 = iup_XkeyAlt
const K_mF7 = iup_XkeyAlt
const K_mF8 = iup_XkeyAlt
const K_mF9 = iup_XkeyAlt
const K_mF10 = iup_XkeyAlt
const K_mF11 = iup_XkeyAlt
const K_mF12 = iup_XkeyAlt
const K_mPrint = iup_XkeyAlt
const K_mMenu = iup_XkeyAlt
const K_yHOME = iup_XkeySys
const K_yUP = iup_XkeySys
const K_yPGUP = iup_XkeySys
const K_yLEFT = iup_XkeySys
const K_yMIDDLE = iup_XkeySys
const K_yRIGHT = iup_XkeySys
const K_yEND = iup_XkeySys
const K_yDOWN = iup_XkeySys
const K_yPGDN = iup_XkeySys
const K_yINS = iup_XkeySys
const K_yDEL = iup_XkeySys
const K_ySP = iup_XkeySys
const K_yTAB = iup_XkeySys
const K_yCR = iup_XkeySys
const K_yBS = iup_XkeySys
const K_yPAUSE = iup_XkeySys
const K_yESC = iup_XkeySys
const K_yCcedilla = iup_XkeySys
const K_yF1 = iup_XkeySys
const K_yF2 = iup_XkeySys
const K_yF3 = iup_XkeySys
const K_yF4 = iup_XkeySys
const K_yF5 = iup_XkeySys
const K_yF6 = iup_XkeySys
const K_yF7 = iup_XkeySys
const K_yF8 = iup_XkeySys
const K_yF9 = iup_XkeySys
const K_yF10 = iup_XkeySys
const K_yF11 = iup_XkeySys
const K_yF12 = iup_XkeySys
const K_yPrint = iup_XkeySys
const K_yMenu = iup_XkeySys
const K_sPlus = iup_XkeyShift
const K_sComma = iup_XkeyShift
const K_sMinus = iup_XkeyShift
const K_sPeriod = iup_XkeyShift
const K_sSlash = iup_XkeyShift
const K_sAsterisk = iup_XkeyShift
const K_cA = iup_XkeyCtrl
const K_cB = iup_XkeyCtrl
const K_cC = iup_XkeyCtrl
const K_cD = iup_XkeyCtrl
const K_cE = iup_XkeyCtrl
const K_cF = iup_XkeyCtrl
const K_cG = iup_XkeyCtrl
const K_cH = iup_XkeyCtrl
const K_cI = iup_XkeyCtrl
const K_cJ = iup_XkeyCtrl
const K_cK = iup_XkeyCtrl
const K_cL = iup_XkeyCtrl
const K_cM = iup_XkeyCtrl
const K_cN = iup_XkeyCtrl
const K_cO = iup_XkeyCtrl
const K_cP = iup_XkeyCtrl
const K_cQ = iup_XkeyCtrl
const K_cR = iup_XkeyCtrl
const K_cS = iup_XkeyCtrl
const K_cT = iup_XkeyCtrl
const K_cU = iup_XkeyCtrl
const K_cV = iup_XkeyCtrl
const K_cW = iup_XkeyCtrl
const K_cX = iup_XkeyCtrl
const K_cY = iup_XkeyCtrl
const K_cZ = iup_XkeyCtrl
const K_c1 = iup_XkeyCtrl
const K_c2 = iup_XkeyCtrl
const K_c3 = iup_XkeyCtrl
const K_c4 = iup_XkeyCtrl
const K_c5 = iup_XkeyCtrl
const K_c6 = iup_XkeyCtrl
const K_c7 = iup_XkeyCtrl
const K_c8 = iup_XkeyCtrl
const K_c9 = iup_XkeyCtrl
const K_c0 = iup_XkeyCtrl
const K_cPlus = iup_XkeyCtrl
const K_cComma = iup_XkeyCtrl
const K_cMinus = iup_XkeyCtrl
const K_cPeriod = iup_XkeyCtrl
const K_cSlash = iup_XkeyCtrl
const K_cSemicolon = iup_XkeyCtrl
const K_cEqual = iup_XkeyCtrl
const K_cBracketleft = iup_XkeyCtrl
const K_cBracketright = iup_XkeyCtrl
const K_cBackslash = iup_XkeyCtrl
const K_cAsterisk = iup_XkeyCtrl
const K_mA = iup_XkeyAlt
const K_mB = iup_XkeyAlt
const K_mC = iup_XkeyAlt
const K_mD = iup_XkeyAlt
const K_mE = iup_XkeyAlt
const K_mF = iup_XkeyAlt
const K_mG = iup_XkeyAlt
const K_mH = iup_XkeyAlt
const K_mI = iup_XkeyAlt
const K_mJ = iup_XkeyAlt
const K_mK = iup_XkeyAlt
const K_mL = iup_XkeyAlt
const K_mM = iup_XkeyAlt
const K_mN = iup_XkeyAlt
const K_mO = iup_XkeyAlt
const K_mP = iup_XkeyAlt
const K_mQ = iup_XkeyAlt
const K_mR = iup_XkeyAlt
const K_mS = iup_XkeyAlt
const K_mT = iup_XkeyAlt
const K_mU = iup_XkeyAlt
const K_mV = iup_XkeyAlt
const K_mW = iup_XkeyAlt
const K_mX = iup_XkeyAlt
const K_mY = iup_XkeyAlt
const K_mZ = iup_XkeyAlt
const K_m1 = iup_XkeyAlt
const K_m2 = iup_XkeyAlt
const K_m3 = iup_XkeyAlt
const K_m4 = iup_XkeyAlt
const K_m5 = iup_XkeyAlt
const K_m6 = iup_XkeyAlt
const K_m7 = iup_XkeyAlt
const K_m8 = iup_XkeyAlt
const K_m9 = iup_XkeyAlt
const K_m0 = iup_XkeyAlt
const K_mPlus = iup_XkeyAlt
const K_mComma = iup_XkeyAlt
const K_mMinus = iup_XkeyAlt
const K_mPeriod = iup_XkeyAlt
const K_mSlash = iup_XkeyAlt
const K_mSemicolon = iup_XkeyAlt
const K_mEqual = iup_XkeyAlt
const K_mBracketleft = iup_XkeyAlt
const K_mBracketright = iup_XkeyAlt
const K_mBackslash = iup_XkeyAlt
const K_mAsterisk = iup_XkeyAlt
const K_yA = iup_XkeySys
const K_yB = iup_XkeySys
const K_yC = iup_XkeySys
const K_yD = iup_XkeySys
const K_yE = iup_XkeySys
const K_yF = iup_XkeySys
const K_yG = iup_XkeySys
const K_yH = iup_XkeySys
const K_yI = iup_XkeySys
const K_yJ = iup_XkeySys
const K_yK = iup_XkeySys
const K_yL = iup_XkeySys
const K_yM = iup_XkeySys
const K_yN = iup_XkeySys
const K_yO = iup_XkeySys
const K_yP = iup_XkeySys
const K_yQ = iup_XkeySys
const K_yR = iup_XkeySys
const K_yS = iup_XkeySys
const K_yT = iup_XkeySys
const K_yU = iup_XkeySys
const K_yV = iup_XkeySys
const K_yW = iup_XkeySys
const K_yX = iup_XkeySys
const K_yY = iup_XkeySys
const K_yZ = iup_XkeySys
const K_y1 = iup_XkeySys
const K_y2 = iup_XkeySys
const K_y3 = iup_XkeySys
const K_y4 = iup_XkeySys
const K_y5 = iup_XkeySys
const K_y6 = iup_XkeySys
const K_y7 = iup_XkeySys
const K_y8 = iup_XkeySys
const K_y9 = iup_XkeySys
const K_y0 = iup_XkeySys
const K_yPlus = iup_XkeySys
const K_yComma = iup_XkeySys
const K_yMinus = iup_XkeySys
const K_yPeriod = iup_XkeySys
const K_ySlash = iup_XkeySys
const K_ySemicolon = iup_XkeySys
const K_yEqual = iup_XkeySys
const K_yBracketleft = iup_XkeySys
const K_yBracketright = iup_XkeySys
const K_yBackslash = iup_XkeySys
const K_yAsterisk = iup_XkeySys

=#


text:
function L_free_particle(mass)
    function L_apply(locals::AbstractLocalTuple)
        v = velocity(locals)
        1//2 * mass * rdot(v,v)
    end
end

function L_harmonic_spring(mass, spring)
    function L_apply(locals::AbstractLocalTuple)
        q = position(locals)
        qv = velocity(locals)
        1//2 * mass * rdot(qv, qv) - 1//2 * spring * rdot(q, q)
    end
end

function L_kepler(particlemass, planetmass)
    function L_apply(locals::AbstractLocalTuple)
        q = position(locals)
        qv = velocity(locals)
        1//2 * particlemass * rdot(qv, qv) + planetmass / sqrt(rdot(q, q))
    end
end

max_stars_repo_name: zoemcc/SICMModelingToolkit

max_issues_repo_name: zoemcc/SICMModelingToolkit

max_forks_repo_name: zoemcc/SICMModelingToolkit


text:
struct Custom
end
JSServe.jsrender(custom::Custom) = DOM.div("i'm a custom struct")

stripnl(x) = strip(x, '\n')
@testset "difflist" begin

    function test_handler(session, request)
        return DOM.div(JSServe.DiffList([Custom(), md"jo", DOM.div("span span span")], dataTestId="difflist"))
    end

    testsession(test_handler) do app
        js_list = query_testid(app, "difflist")
        difflist = children(app.dom)[1]

        connected_list = Observable(Any[])
        JSServe.connect!(connected_list, difflist)
        test_connection() = connected_list[] == difflist.children

        @wait_for evaljs(app, js_list.children.length) == 3
        @test evaljs(app, js"$(js_list).children[0].innerText") == "i'm a custom struct"
        @test stripnl(evaljs(app, js"$(js_list).children[1].innerText")) == "jo"
        @test evaljs(app, js"$(js_list).children[2].innerText") == "span span span"

        empty!(difflist)
        @wait_for evaljs(app, js_list.children.length) == 0
        push!(difflist, md"new node 1")
        push!(difflist, md"new node 2")
        push!(difflist, md"new node 3")

        @test test_connection()

        @test stripnl(evaljs(app, js"$(js_list).children[0].innerText")) == "new node 1"
        @test stripnl(evaljs(app, js"$(js_list).children[1].innerText")) == "new node 2"
        @test stripnl(evaljs(app, js"$(js_list).children[2].innerText")) == "new node 3"
        append!(difflist, [md"append", md"much fun"])
        @test stripnl(evaljs(app, js"$(js_list).children[3].innerText")) == "append"
        @test stripnl(evaljs(app, js"$(js_list).children[4].innerText")) == "much fun"

        difflist[2] = md"old node 2"
        @test stripnl(evaljs(app, js"$(js_list).children[1].innerText")) == "old node 2"

        JSServe.replace_children(difflist, [md"meowdy", md"monsieur", md"maunz"])

        @wait_for evaljs(app, js_list.children.length) == 3

        @test test_connection()

        @test stripnl(evaljs(app, js"$(js_list).children[0].innerText")) == "meowdy"
        @test stripnl(evaljs(app, js"$(js_list).children[1].innerText")) == "monsieur"
        @test stripnl(evaljs(app, js"$(js_list).children[2].innerText")) == "maunz"

        insert!(difflist, 2, md"rowdy")
        @test stripnl(evaljs(app, js"$(js_list).children[1].innerText")) == "rowdy"
        delete!(difflist, [1, 4])

        @test test_connection()

        @test stripnl(evaljs(app, js"$(js_list).children[0].innerText")) == "rowdy"
        @test stripnl(evaljs(app, js"$(js_list).children[1].innerText")) == "monsieur"
        @test children(difflist) == [md"rowdy", md"monsieur"]

        @test test_connection()
        JSServe.replace_children(difflist, [md"$i" for i in 1:3000])
        insert!(difflist, 2999, md"hehe")
        @test length(difflist) == 3001
        @test difflist[2999] == md"hehe"
        @test difflist[3000].content[1] == 2999
        @test difflist[3001].content[1] == 3000
        @test stripnl(evaljs(app, js"$(js_list).children[2999-1].innerText")) == "hehe"
        @test stripnl(evaljs(app, js"$(js_list).children[3000-1].innerText")) == "2999"
        @test stripnl(evaljs(app, js"$(js_list).children[3001-1].innerText")) == "3000"
    end
end


text:
# ---
# title: 419. Battleships in a Board
# id: problem419
# author: Tian Jun
# date: 2020-10-31
# difficulty: Medium
# categories:
# link: <https://leetcode.com/problems/battleships-in-a-board/description/>
# hidden: true
# ---
#
# Given an 2D board, count how many battleships are in it. The battleships are
# represented with `'X'`s, empty slots are represented with `'.'`s. You may
# assume the following rules:
#
#   * You receive a valid board, made of only battleships or empty slots.
#   * Battleships can only be placed horizontally or vertically. In other words, they can only be made of the shape `1xN` (1 row, N columns) or `Nx1` (N rows, 1 column), where N can be of any size.
#   * At least one horizontal or vertical cell separates between two battleships - there are no adjacent battleships.
#
# **Example:**
#
#
#
#     X..X
#     ...X
#     ...X
#
#
# In the above board there are 2 battleships.
#
# **Invalid Example:**
#
#
#
#     ...X
#     XXXX
#     ...X
#
#
# This is an invalid board that you will not receive - as battleships will
# always have a cell separating between them.
#
# **Follow up:**
# Could you do it in **one-pass** , using only **O(1) extra memory** and
# **without modifying** the value of the board?
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end


text:
module BatchedBLAS

using CUDA

export batched_dot!
export batched_gemv!, batched_symv!, batched_spmv!
export batched_ger!, batched_syr!, batched_spr!

const IntOrFloat = Union{Integer,AbstractFloat}

# https://github.com/JuliaLang/julia/issues/40469
maybe_cast(::Type, x) = x
maybe_cast(::Type{T}, x::AbstractFloat) where T<:Integer =
        round(T, clamp(x, typemin(T), typemax(T)))

function configurator(config, dims)
    if length(dims)==1
        xthreads = min(32, dims[1])
        xblocks = cld(dims[1], xthreads)
        return (xthreads,), (xblocks,)
    elseif length(dims)==2
        xthreads = min(32, dims[1])
        ythreads = min(fld(config.threads, xthreads), cld(prod(dims), xthreads))
        xblocks = cld(dims[1], xthreads)
        yblocks = cld(dims[2], ythreads)
        return (xthreads, ythreads), (xblocks, yblocks)
    end
end

"""
    batched_dot!(o, x, y)

In-place batched vector-vector multiplication, equivalent to `o[:,k] =
transpose(x[:,k]) * y[:,k]` for all `k`.  All inputs can have eltypes of
either AbstractFloats or Integers.
"""
function batched_dot!(o::CuVector{To}, x::CuMatrix{Tx}, y::CuMatrix{Ty}) where {
                          To<:IntOrFloat, Tx<:IntOrFloat, Ty<:IntOrFloat}

    function kernel(T, o, x, y)
        k = threadIdx().x + (blockIdx().x - 1) * blockDim().x

        @inbounds if k<=size(x,2)
            o[k] = 0
            for i=1:size(x,1)
                o[k] += maybe_cast(To, x[i,k] * y[i,k])
            end
        end
    end

    T = promote_type(To, Tx, Ty)
    kernel = @cuda name="batched_dot!" launch=false kernel(T, o, x, y)
    config = launch_configuration(kernel.fun)
    threads, blocks = configurator(config, (size(o,1),))
    kernel(T, o, x, y; threads=threads, blocks=blocks)
end

"""
    batched_gemv!(tA, alpha, A, x, beta, y)

In-place batched matrix-vector multiplication and addition, equivalent
to `y[:,k] = alpha[k] * A[:,:,k] * x[:,k] + beta[k] * y[:,k]` for all
`k`. `A` can optionally be transposed with `tA` as `N`, `T`, or `C`.
All other inputs can have eltypes of either AbstractFloats or Integers.
`alpha` and `beta` can also be scalars.
"""
function batched_gemv!(tA::AbstractChar,
                       alpha::Talpha, A::CuArray{TA,3}, x::CuMatrix{Tx},
                       beta::Tbeta, y::CuMatrix{Ty}) where {
                           Talpha<:Union{IntOrFloat, CuVector{<:IntOrFloat}},
                           TA<:IntOrFloat, Tx<:IntOrFloat,
                           Tbeta<:Union{IntOrFloat, CuVector{<:IntOrFloat}},
                           Ty<:IntOrFloat}

    function kernel(T, tA, alpha, A, x, beta, y)
        i = threadIdx().x + (blockIdx().x - 1) * blockDim().x
        k = threadIdx().y + (blockIdx().y - 1) * blockDim().y

        if tA=='N'
            @inbounds if k<=size(y,2) && i<=size(y,1)
                tmp = T(0)
                for j=1:size(x,1)
                    tmp += A[i,j,k] * x[j,k]
                end
                thisalpha = Talpha<:CuVector ? alpha[k] : alpha
                thisbeta = Tbeta<:CuVector ? beta[k] : beta
                y[i,k] = maybe_cast(Ty, thisalpha*tmp + thisbeta*y[i,k])
            end
        elseif tA=='T'
            @inbounds if k<=size(y,2) && i<=size(y,1)
                tmp = T(0)
                for j=1:size(x,1)
                    tmp += A[j,i,k] * x[j,k]
                end
                thisalpha = Talpha<:CuVector ? alpha[k] : alpha
                thisbeta = Tbeta<:CuVector ? beta[k] : beta
                y[i,k] = maybe_cast(Ty, thisalpha*tmp + thisbeta*y[i,k])
            end
        elseif tA=='C'
            @inbounds if k<=size(y,2) && i<=size(y,1)
                tmp = T(0)
                for j=1:size(x,1)
                    tmp += adjoint(A[j,i,k]) * x[j,k]
                end
                thisalpha = Talpha<:CuVector ? alpha[k] : alpha
                thisbeta = Tbeta<:CuVector ? beta[k] : beta
                y[i,k] = maybe_cast(Ty, thisalpha*tmp + thisbeta*y[i,k])
            end
        else
            throw(ArgumentError("`tA` should be 'N', 'T', or 'C'"))
        end
        return nothing
    end

    T = promote_type(TA, Tx, Ty)
    kernel = @cuda name="batched_gemv!" launch=false kernel(T, tA, alpha, A, x, beta, y)
    config = launch_configuration(kernel.fun)
    threads, blocks = configurator(config, (size(y,1),size(y,2)))
    kernel(T, tA, alpha, A, x, beta, y; threads=threads, blocks=blocks)
end

"""
    batched_symv!(uplo, alpha, A, x, beta, y)

In-place batched matrix-vector multiplication and addition, equivalent to
`y[:,k] = alpha[k] * A[:,:,k] * x[:,k] + beta[k] * y[:,k]` for all `k`.  `A`
is assumed to be symmetric.  Only the `uplo` (either 'U' or 'L') triangle of
`A` is used.  All other inputs can have eltypes of either AbstractFloats
or Integers.  `alpha` and `beta` can also be scalars.
"""
function batched_symv!(uplo::AbstractChar,
                       alpha::Talpha, A::CuArray{TA,3}, x::CuMatrix{Tx},
                       beta::Tbeta, y::CuMatrix{Ty}) where {
                           Talpha<:Union{IntOrFloat, CuVector{<:IntOrFloat}},
                           TA<:IntOrFloat, Tx<:IntOrFloat,
                           Tbeta<:Union{IntOrFloat, CuVector{<:IntOrFloat}},
                           Ty<:IntOrFloat}

    function kernel(T, uplo, alpha, A, x, beta, y)
        i = threadIdx().x + (blockIdx().x - 1) * blockDim().x
        k = threadIdx().y + (blockIdx().y - 1) * blockDim().y

        if uplo=='U'
            @inbounds if k<=size(y,2) && i<=size(y,1)
                tmp = T(0)
                for j=1:size(x,1)
                    ijmin,ijmax = minmax(i,j)
                    tmp += A[ijmin,ijmax,k] * x[j,k]
                end
                thisalpha = Talpha<:CuVector ? alpha[k] : alpha
                thisbeta = Tbeta<:CuVector ? beta[k] : beta
                y[i,k] = maybe_cast(Ty, thisalpha*tmp + thisbeta*y[i,k])
            end
        elseif uplo=='L'
            @inbounds if k<=size(y,2) && i<=size(y,1)
                tmp = T(0)
                for j=1:size(x,1)
                    ijmin,ijmax = minmax(i,j)
                    tmp += A[ijmax,ijmin,k] * x[j,k]
                end
                thisalpha = Talpha<:CuVector ? alpha[k] : alpha
                thisbeta = Tbeta<:CuVector ? beta[k] : beta
                y[i,k] = maybe_cast(Ty, thisalpha*tmp + thisbeta*y[i,k])
            end
        else
            throw(ArgumentError("`uplo` should be 'U' or 'L'"))
        end
        return nothing
    end

    T = promote_type(TA, Tx, Ty)
    kernel = @cuda name="batched_symv!" launch=false kernel(T, uplo, alpha, A, x, beta, y)
    config = launch_configuration(kernel.fun)
    threads, blocks = configurator(config, (size(y,1),size(y,2)))
    kernel(T, uplo, alpha, A, x, beta, y; threads=threads, blocks=blocks)
end

"""
  batched_spmv!(ul, alpha, A, x, beta, y)

In-place batched matrix-vector multiplication and addition, equivalent to
`y[:,k] = alpha[k] * A[:,:,k] * x[:,k] + beta[k] * y[:,k]` for all `k`.
`A` must be symmetric, and `uplo` specifies whether the upper ('U') or lower
('L') triangle was packed.  All other inputs can have eltypes of either
AbstractFloats or Integers.  `alpha` and `beta` can also be scalars.
"""
function batched_spmv!(uplo::AbstractChar,
                       alpha::Talpha, A::CuMatrix{TA}, x::CuMatrix{Tx},
                       beta::Tbeta, y::CuMatrix{Ty}) where {
                         Talpha<:Union{IntOrFloat, CuVector{<:IntOrFloat}},
                         TA<:IntOrFloat, Tx<:IntOrFloat,
                         Tbeta<:Union{IntOrFloat, CuVector{<:IntOrFloat}},
                         Ty<:IntOrFloat}

    function kernel(T, uplo, alpha, A, x, beta, y)
        i = threadIdx().x + (blockIdx().x - 1) * blockDim().x
        k = threadIdx().y + (blockIdx().y - 1) * blockDim().y

        if uplo=='U'
            @inbounds if k<=size(y,2) && i<=size(y,1)
                tmp = T(0)
                for j=1:size(x,1)
                    ijmin,ijmax = minmax(i,j)
                    h = ijmin+(ijmax*(ijmax-1))>>1
                    tmp += A[h,k] * x[j,k]
                end
                thisalpha = Talpha<:CuVector ? alpha[k] : alpha
                thisbeta = Tbeta<:CuVector ? beta[k] : beta
                y[i,k] = maybe_cast(Ty, thisalpha*tmp + thisbeta*y[i,k])
            end
        elseif uplo=='L'
            n = round(Int, (sqrt(8*size(A,1))-1)/2)
            @inbounds if k<=size(y,2) && i<=size(y,1)
                tmp = T(0)
                for j=1:size(x,1)
                    ijmin,ijmax = minmax(i,j)
                    h = ijmax+((2n-ijmin)*(ijmin-1))>>1
                    tmp += A[h,k] * x[j,k]
                end
                thisalpha = Talpha<:CuVector ? alpha[k] : alpha
                thisbeta = Tbeta<:CuVector ? beta[k] : beta
                y[i,k] = maybe_cast(Ty, thisalpha*tmp + thisbeta*y[i,k])
            end
        else
            throw(ArgumentError("`uplo` should be 'U' or 'L'"))
        end
        return nothing
    end

    T = promote_type(TA, Tx, Ty)
    kernel = @cuda name="batched_spmv!" launch=false kernel(T, uplo, alpha, A, x, beta, y)
    config = launch_configuration(kernel.fun)
    threads, blocks = configurator(config, (size(y,1),size(y,2)))
    kernel(T, uplo, alpha, A, x, beta, y; threads=threads, blocks=blocks)
end

"""
    batched_ger!(alpha, x, y, A)

In-place rank-1 update of matrix `A` with vectors `x` and `y` as `alpha[k] *
x[:,k] * transpose(y[:,k]) + A[:,:,k]` for all `k`.  All inputs can have
eltypes of either AbstractFloats or Integers.  `alpha` can also be a scalar.
"""
function batched_ger!(alpha::Talpha, x::CuMatrix{Tx}, y::CuMatrix{Ty}, A::CuArray{TA,3}) where {
                          Talpha<:Union{IntOrFloat, CuVector{<:IntOrFloat}},
                          Tx<:IntOrFloat, Ty<:IntOrFloat, TA<:IntOrFloat}

    function kernel(alpha, x, y, A)
        i = threadIdx().x + (blockIdx().x - 1) * blockDim().x
        k = threadIdx().y + (blockIdx().y - 1) * blockDim().y

        @inbounds if k<=size(x,2) && i<=size(x,1)
            for j=1:size(x,1)
                thisalpha = Talpha<:CuVector ? alpha[k] : alpha
                A[i,j,k] += maybe_cast(TA, thisalpha * x[i,k] * y[j,k])
            end
        end
    end

    kernel = @cuda name="batched_ger_vector!" launch=false kernel(alpha, x, y, A)
    config = launch_configuration(kernel.fun)
    threads, blocks = configurator(config, (size(x,1),size(x,2)))
    kernel(alpha, x, y, A; threads=threads, blocks=blocks)
end

"""
    batched_syr!(uplo, alpha, x, A)

In-place rank-1 update of symmetric matrix `A` with vector `x` as `alpha[k]
* x[:,k] * transpose(x[:,k]) + A[:,:,k]` for all `k`.  `A` is assumed to
be symmetric.  Only the `uplo` (either 'U' or 'L') triangle of `A` is used.
All other inputs can have eltypes of either AbstractFloats or Integers.
`alpha` can also be a scalar.
"""
function batched_syr!(uplo::AbstractChar,
                      alpha::Talpha, x::CuMatrix{Tx}, A::CuArray{TA,3}) where {
                          Talpha<:Union{IntOrFloat, CuVector{<:IntOrFloat}},
                          Tx<:IntOrFloat, TA<:IntOrFloat}

    function kernel(uplo, alpha, x, A)
        i = threadIdx().x + (blockIdx().x - 1) * blockDim().x
        k = threadIdx().y + (blockIdx().y - 1) * blockDim().y

        if uplo=='U'
            @inbounds if k<=size(x,2) && i<=size(x,1)
                for j=size(x,1):-1:1
                    j<i && break
                    thisalpha = Talpha<:CuVector ? alpha[k] : alpha
                    A[i,j,k] += maybe_cast(TA, thisalpha * x[i,k] * x[j,k])
                end
            end
        elseif uplo=='L'
            @inbounds if k<=size(x,2) && i<=size(x,1)
                for j=1:size(x,1)
                    j>i && break
                    thisalpha = Talpha<:CuVector ? alpha[k] : alpha
                    A[i,j,k] += maybe_cast(TA, thisalpha * x[i,k] * x[j,k])
                end
            end
        else
            throw(ArgumentError("`uplo` should be 'U' or 'L'"))
        end
    end

    kernel = @cuda name="batched_syr_vector!" launch=false kernel(uplo, alpha, x, A)
    config = launch_configuration(kernel.fun)
    threads, blocks = configurator(config, (size(x,1),size(x,2)))
    kernel(uplo, alpha, x, A; threads=threads, blocks=blocks)
end

"""
    batched_spr!(uplo, alpha, x, A)

In-place rank-1 update of packed symmetric matrix `A` with vector `x`
as `alpha[k] * x[:,k] * transpose(x[:,k]) + A[:,:,k]` for all `k`.  `A`
must be symmetric, and `uplo` specifies whether the upper ('U') or lower
('L') triangle was packed.  All other inputs can have eltypes of either
AbstractFloats or Integers.  `alpha` can also be a scalar.
"""
function batched_spr!(uplo::AbstractChar,
                      alpha::Talpha, x::CuMatrix{Tx}, A::CuMatrix{TA}) where {
                          Talpha<:Union{IntOrFloat, CuVector{<:IntOrFloat}},
                          Tx<:IntOrFloat, TA<:IntOrFloat}

    function kernel(uplo, alpha, x, A)
        i = threadIdx().x + (blockIdx().x - 1) * blockDim().x
        k = threadIdx().y + (blockIdx().y - 1) * blockDim().y

        if uplo=='U'
            @inbounds if k<=size(x,2) && i<=size(x,1)
                for j=size(x,1):-1:1
                    j<i && break
                    h = i+(j*(j-1))>>1
                    thisalpha = Talpha<:CuVector ? alpha[k] : alpha
                    A[h,k] += maybe_cast(TA, thisalpha * x[i,k] * x[j,k])
                end
            end
        elseif uplo=='L'
            n = round(Int, (sqrt(8*size(A,1))-1)/2)
            @inbounds if k<=size(x,2) && i<=size(x,1)
                for j=1:size(x,1)
                    j>i && break
                    h = i+((2n-j)*(j-1))>>1
                    thisalpha = Talpha<:CuVector ? alpha[k] : alpha
                    A[h,k] += maybe_cast(TA, thisalpha * x[i,k] * x[j,k])
                end
            end
        else
            throw(ArgumentError("`uplo` should be 'U' or 'L'"))
        end
    end

    kernel = @cuda name="batched_spr_vector!" launch=false kernel(uplo, alpha, x, A)
    config = launch_configuration(kernel.fun)
    threads, blocks = configurator(config, (size(x,1),size(x,2)))
    kernel(uplo, alpha, x, A; threads=threads, blocks=blocks)
end

end


text:
export BrickletDistanceUSDistanceCallbackThreshold
struct BrickletDistanceUSDistanceCallbackThreshold
    option::Char
    min::Integer
    max::Integer
end

export BrickletDistanceUSIdentity
struct BrickletDistanceUSIdentity
    uid::String
    connected_uid::String
    position::Char
    hardware_version::Vector{Integer}
    firmware_version::Vector{Integer}
    device_identifier::Integer
end

export BrickletDistanceUS
"""
Measures distance between 2cm and 400cm with ultrasound
"""
mutable struct BrickletDistanceUS <: TinkerforgeDevice
    ipcon::IPConnection
    deviceInternal::PyObject

    """
    Creates an object with the unique device ID *uid* and adds it to
    the IP Connection *ipcon*.
    """
    function BrickletDistanceUS(uid::String, ipcon::IPConnection)
        package = pyimport("tinkerforge.bricklet_distance_us")
        deviceInternal = package.BrickletDistanceUS(uid, ipcon.ipconInternal)

        return new(ipcon, deviceInternal)
    end
end

export get_distance_value
"""
    $(SIGNATURES)

Returns the current distance value measured by the sensor.
A small value corresponds to a small distance, a big
value corresponds to a big distance. The relation between the measured distance
value and the actual distance is affected by the 5V supply voltage (deviations
in the supply voltage result in deviations in the distance values) and is
non-linear (resolution is bigger at close range).

If you want to get the distance value periodically, it is recommended to
use the :cb:`Distance` callback and set the period with
:func:`Set Distance Callback Period`.
"""
function get_distance_value(device::BrickletDistanceUS)
    return device.deviceInternal.get_distance_value()
end

export set_distance_callback_period
"""
    $(SIGNATURES)

Sets the period with which the :cb:`Distance` callback is triggered
periodically. A value of 0 turns the callback off.

Der :cb:`Distance` callback is only triggered if the distance value has changed
since the last triggering.
"""
function set_distance_callback_period(device::BrickletDistanceUS, period)
    device.deviceInternal.set_distance_callback_period(period)
end

export get_distance_callback_period
"""
    $(SIGNATURES)

Returns the period as set by :func:`Set Distance Callback Period`.
"""
function get_distance_callback_period(device::BrickletDistanceUS)
    return device.deviceInternal.get_distance_callback_period()
end

export set_distance_callback_threshold
"""
    $(SIGNATURES)

Sets the thresholds for the :cb:`Distance Reached` callback.

The following options are possible:

.. csv-table::
 :header: "Option", "Description"
 :widths: 10, 100

 "'x'",    "Callback is turned off"
 "'o'",    "Callback is triggered when the distance value is *outside* the min and max values"
 "'i'",    "Callback is triggered when the distance value is *inside* the min and max values"
 "'<'",    "Callback is triggered when the distance value is smaller than the min value (max is ignored)"
 "'>'",    "Callback is triggered when the distance value is greater than the min value (max is ignored)"
"""
function set_distance_callback_threshold(device::BrickletDistanceUS, option, min, max)
    device.deviceInternal.set_distance_callback_threshold(option, min, max)
end

export get_distance_callback_threshold
"""
    $(SIGNATURES)

Returns the threshold as set by :func:`Set Distance Callback Threshold`.
"""
function get_distance_callback_threshold(device::BrickletDistanceUS)
    return device.deviceInternal.get_distance_callback_threshold()
end

export set_debounce_period
"""
    $(SIGNATURES)

Sets the period with which the threshold callbacks

* :cb:`Distance Reached`,

are triggered, if the thresholds

* :func:`Set Distance Callback Threshold`,

keep being reached.
"""
function set_debounce_period(device::BrickletDistanceUS, debounce)
    device.deviceInternal.set_debounce_period(debounce)
end

export get_debounce_period
"""
    $(SIGNATURES)

Returns the debounce period as set by :func:`Set Debounce Period`.
"""
function get_debounce_period(device::BrickletDistanceUS)
    return device.deviceInternal.get_debounce_period()
end

export set_moving_average
"""
    $(SIGNATURES)

Sets the length of a `moving averaging <https://en.wikipedia.org/wiki/Moving_average>`__
for the distance value.

Setting the length to 0 will turn the averaging completely off. With less
averaging, there is more noise on the data.
"""
function set_moving_average(device::BrickletDistanceUS, average)
    device.deviceInternal.set_moving_average(average)
end

export get_moving_average
"""
    $(SIGNATURES)

Returns the length moving average as set by :func:`Set Moving Average`.
"""
function get_moving_average(device::BrickletDistanceUS)
    return device.deviceInternal.get_moving_average()
end

export get_identity
"""
    $(SIGNATURES)

Returns the UID, the UID where the Bricklet is connected to,
the position, the hardware and firmware version as well as the
device identifier.

The position can be 'a', 'b', 'c', 'd', 'e', 'f', 'g' or 'h' (Bricklet Port).
A Bricklet connected to an :ref:`Isolator Bricklet <isolator_bricklet>` is always at
position 'z'.

The device identifier numbers can be found :ref:`here <device_identifier>`.
|device_identifier_constant|
"""
function get_identity(device::BrickletDistanceUS)
    return device.deviceInternal.get_identity()
end

export register_callback
"""
Registers the given *function* with the given *callback_id*.
"""
function register_callback(device::BrickletDistanceUS, callback_id, function_)
    if isnothing(function_)
        device.registered_callbacks.pop(callback_id, None)
    else
        device.registered_callbacks[callback_id] = function_
    end
end


text:
"""
Script for comparison of two versions of  SdErgmDirBin0Rec0, i.e. one parameter for total number of links and one for reciprocity, a version based on MLE and one on PMLE

To Do:

"""

using Statistics
using PyPlot
pygui(true)

using ScoreDrivenERGM
using ScoreDrivenERGM.Utilities
using ScoreDrivenERGM.Scalings

using StaticNets:ergm_par_from_mean_vals
using DynNets:fooSdErgmDirBin0Rec0_mle, sample_est_mle_pmle, dgp_missp, SdErgmDirBin0Rec0_mle

import StaticNets: ErgmDirBin0Rec0(), ergm_par_from_mean_vals,exp_val_stats, diadProbFromPars, samplSingMatCan, stats_from_mat

model_mle = fooSdErgmDirBin0Rec0_mle
staModel = ErgmDirBin0Rec0()

## Search an appropriate dgp, i.e. one that stays away form physical bounds

#The following could be updated using sample_time_var_par_from_dgp()
function scaled_sin_sample_est_mle_pmle_var_N(alphaMeanVal::Function, N, T, nSample; regimeString ="", plotDgpOrigin=false, plotFlag=false)

    betaMeanVal(N) = alphaMeanVal(N)/5

    amplAlpha = 0.3
    amplBeta = 0.5

    nCycles = 2

    alpha_beta_seq_given_N(N) = hcat(dgpSin( (1-amplAlpha/2)*alphaMeanVal(N), (1+amplAlpha/2)*alphaMeanVal(N), nCycles, T; phase = 10*rand()), dgpSin( (1-amplBeta/2)*betaMeanVal(N), (1+amplBeta/2)*betaMeanVal(N), nCycles, T; phase = 10*rand()))

    alpha_beta_to_theta_eta(α, β, N) = collect(ergm_par_from_mean_vals(staModel, α*n_pox_dir_links(N), β*n_pox_dir_links(N), N))

    theta_eta_to_alpha_beta(θ, η, N) =  collect(exp_val_stats(staModel, θ, η, N))./n_pox_dir_links(N)

    get_theta_eta_seq_from_alpha_beta(N) = reduce(hcat, [alpha_beta_to_theta_eta(ab[1], ab[2], N) for ab in eachrow(alpha_beta_seq_given_N(N))])

    if plotDgpOrigin
        startAlphaBetaSeq = alpha_beta_seq_given_N(N)
        thetaEtaSeq = permutedims(reduce(hcat, [alpha_beta_to_theta_eta(ab[1], ab[2], N) for ab in eachrow(startAlphaBetaSeq)]))

        fig, ax = subplots(2,2)
        ax[1,1].plot(startAlphaBetaSeq[:,1])
        ax[1,1].set_title("Alpha Dgp")
        ax[2,1].plot(startAlphaBetaSeq[:,2])
        ax[2,1].set_title("Beta Dgp")
        ax[1,2].plot(thetaEtaSeq[:,1])
        ax[1,2].set_title("Theta Dgp")
        ax[2,2].plot(thetaEtaSeq[:,2])
        ax[2,2].set_title("Eta Dgp")
        fig.tight_layout()
        fig.suptitle("$regimeString Dgp for N = $N ")
    end

    # In case we want to add a check before filtering, start from this commented lines
    # minNPairs = 10
    # lBoundBeta = minNPairs/n_pox_dir_links(N)
    # uBoundBeta = alphaMeanVal(N)/2 -  minNPairs/n_pox_dir_links(N)
    # startLRSeq = alpha_beta_seq_given_N(N).*n_pox_dir_links(N)
    # endAlphaBetaSeq = permutedims( reduce(hcat, [theta_eta_to_alpha_beta(θ, η, N) for (θ, η) in eachcol(get_theta_eta_seq_from_alpha_beta(N))]) )
    # θ_0, η_0 = get_theta_eta_seq_from_alpha_beta(N)[:,1]
    # A_vec = [stats_from_mat(staModel, samplSingMatCan(staModel, diadProbFromPars(staModel, [θ_0, η_0]), N)) for i=1:100]

    return sample_est_mle_pmle(model_mle, get_theta_eta_seq_from_alpha_beta(N), N, nSample; plotFlag=plotFlag, regimeString=regimeString)
end

# Check that misspecified filters run without evident issues on few samples at N extrema

denseAlphaScal(N) = denseLScal_DirBin(N)/n_pox_dir_links(N)
sparseAlphaScal(N) = sparseLScal_DirBin(N)/n_pox_dir_links(N)
T = 200
scaled_sin_sample_est_mle_pmle(denseAlphaScal, 10, T, 2, regimeString = "Dense", plotDgpOrigin=true)
scaled_sin_sample_est_mle_pmle(sparseAlphaScal, 10, T, 2, regimeString = "Sparse", plotDgpOrigin=true)

scaled_sin_sample_est_mle_pmle(denseAlphaScal, 100, T, 2,regimeString = "Dense", plotDgpOrigin=true)
scaled_sin_sample_est_mle_pmle(sparseAlphaScal, 100, T, 2, regimeString = "Sparse", plotDgpOrigin=true)

"""
for a single misspecified dgp compare mle and mple rmse for varying network sizes
"""
function scaling_comparison(model_mle::SdErgmDirBin0Rec0_mle, Nvals, alphaScal, nSample, T; plotFlag = false, regimeString="")

    staModel = ErgmDirBin0Rec0()

    rmse_mle = []
    rmse_pmle = []
    mean_vals = []

    for N in Nvals

        res = scaled_sin_sample_est_mle_pmle(alphaScal, N, T, nSample, regimeString=regimeString; plotFlag=plotFlag)

        push!(rmse_pmle, res.rmse_pmle)
        push!(rmse_mle, res.rmse_mle)

    end

    return rmse_mle, rmse_pmle, mean_vals
end

using Serialization

T = 200
# Nvals= [20, 100]# round.(Int, collect(20:20:100))
Nvals=  round.(Int, collect(20:10:100))
nSample = 50

savePath = (@__DIR__ )* "\\revision_data"

scaling_res = scaling_comparison(model_mle, Nvals, sparseAlphaScal, nSample, T; plotFlag = false, regimeString = "Sparse")
serialize(savePath * "\\sparse_nSample$(nSample)_precise_filter_init.jls", scaling_res)

scaling_res =scaling_comparison(model_mle, Nvals, denseAlphaScal, nSample, T; plotFlag = false, regimeString = "Dense")
serialize(savePath * "\\dense_nSample$(nSample)_precise_filter_init.jls", scaling_res)

nSample = 50
rmse_mle_sparse, rmse_pmle_sparse, mean_vals_sparse = deserialize(savePath * "\\sparse_nSample$(nSample)_precise_filter_init.jls")

rmse_mle_dense, rmse_pmle_dense, mean_vals_dense = deserialize(savePath * "\\dense_nSample$(nSample)_precise_filter_init.jls")

fig, ax = subplots(2,2)
ax[1,1].plot(Nvals, reduce(hcat,rmse_mle_sparse)[1,:], "-.b")
ax[1,1].plot(Nvals, reduce(hcat,rmse_pmle_sparse)[1,:], "-.r")
ax[2,1].plot(Nvals, reduce(hcat,rmse_mle_sparse)[2,:], "-.b")
ax[2,1].plot(Nvals, reduce(hcat,rmse_pmle_sparse)[2,:], "-.r")
ax[1,1].set_title("Sparse regime")

ax[1,2].plot(Nvals, reduce(hcat,rmse_mle_dense)[1,:], "-.b")
ax[1,2].plot(Nvals, reduce(hcat,rmse_pmle_dense)[1,:], "-.r")
ax[2,2].plot(Nvals, reduce(hcat,rmse_mle_dense)[2,:], "-.b")
ax[2,2].plot(Nvals, reduce(hcat,rmse_pmle_dense)[2,:], "-.r")
ax[1,2].set_title("Dense  regime")

# N=80
# linkScal(N) = semiDenseLScal_DirBin(N)
# recScal(N) = min(5 + 10 * erdosRenyRecScal_DirBin(linkScal(N), N), linkScal(N)/2 -1)
# percAmpl = 0.5
# θ_0_min, η_0_min = ergm_par_from_mean_vals(staModel, linkScal(N)*(1-percAmpl), recScal(N)*(1-percAmpl), N)
# θ_0_max, η_0_max = ergm_par_from_mean_vals(staModel, linkScal(N)*(1+percAmpl), recScal(N)*(1+percAmpl), N)
# (linkScal(N), recScal(N))

# η_0_minMax = sort([η_0_min, η_0_max])
# θ_0_minMax = sort([θ_0_min, θ_0_max])
# resSin = sample_est_mle_pmle(model_mle, dgp_missp(model_mle, T, θ_0_minMax, η_0_minMax, "sin"), N, nSample)
# resAR = sample_est_mle_pmle(model_mle, dgp_missp(model_mle, T, θ_0_minMax, η_0_minMax, "AR"), N, nSample)

max_stars_repo_name: domenicodigangi/SDERGM

max_issues_repo_name: domenicodigangi/SDERGM

max_forks_repo_name: domenicodigangi/SDERGM


text:
# debug.jl
# this file contain debugging utilities

# global variable that determines debug level
# must be a literal value
global const DB_LEVEL = 1
# global const DB_LEVEL = 1

@doc """
### PDESolver.debug1

  This macro returns the expression if the debug level is greater than or
  equal to 1.

  Inputs/Outputs:
    expr1: an expression to be executed only if the debug level is greater
           than 1

"""->
macro debug1(expr1)
#  println("entered macro do_db")
  if DB_LEVEL >= 1
#    println("at compile time, in DB_Level >=", DB_LEVEL)
    return quote
#      println("runtime expression")
        $(esc(expr1))
    end
  else
    return nothing
  end

end  # end macro

@doc """
### PDESolver.debug2

  This macro returns the expression if the debug level is greater than or
  equal to 2

  Inputs/Outputs:
    expr1: an expression to be executed only if the debug level is greater
           than 2

"""->
macro debug2(expr1)
#  println("entered macro do_db")
  if DB_LEVEL >= 2
#    println("at compile time, in DB_Level >=", DB_LEVEL)
    return quote
#      println("runtime expression")
        $(esc(expr1))
    end
  else
    return nothing
  end

end  # end macro


text:
ext(::MIME"image/png") = ".png"
ext(::MIME"image/svg+xml") = ".svg"

function Base.show(io::IO, m::MIME"image/png", plt::RObject)
  @show f = tempname() * ext(m)
  @show  width, height = get(io, :juno_plotsize,  [6*72, 5*72])

  R"""
    png($f, width=$width, height=$height, units="px")
    show($plt)
    dev.off()
  """
  show(io, m, load(f))
end

function Base.show(io::IO, m::MIME"image/svg+xml", plt::RObject)
  @show f = tempname() * ext(m)
  #@show  width, height = get(io, :juno_plotsize,  [6, 5])

  R"""
    svg($f)
    show($plt)
    dev.off()
  """
  ## Copied from RCall.jl/src/ijulia.jl
  open(f) do f
      r = randstring()
      d = read(f, String)
      d = replace(d, "id=\"glyph" => "id=\"glyph"*r)
      d = replace(d, "href=\"#glyph" => "href=\"#glyph"*r)
      display(m,d)
  end
end

isggplot(plt::RObject) = rcopy(R"""ggplot2::is.ggplot($plt)""")

Base.showable(m::Union{MIME"image/png",MIME"image/svg+xml"}, plt::RObject) = is_ggplot(plt)


text:
# ---
# title: 236. Lowest Common Ancestor of a Binary Tree
# id: problem236
# author: zhwang
# date: 2022-02-10
# difficulty: Medium
# categories: Tree
# link: <https://leetcode.com/problems/lowest-common-ancestor-of-a-binary-tree/description/>
# hidden: true
# ---
#
# Given a binary tree, find the lowest common ancestor (LCA) of two given nodes
# in the tree.
#
# According to the [definition of LCA on
# Wikipedia](https://en.wikipedia.org/wiki/Lowest_common_ancestor): "The lowest
# common ancestor is defined between two nodes p and q as the lowest node in T
# that has both p and q as descendants (where we allow **a node to be a
# descendant of itself** )."
#
#
#
# **Example 1:**
#
# ![](https://assets.leetcode.com/uploads/2018/12/14/binarytree.png)
#
#
#
#     Input: root = [3,5,1,6,2,0,8,null,null,7,4], p = 5, q = 1
#     Output: 3
#     Explanation: The LCA of nodes 5 and 1 is 3.
#
#
# **Example 2:**
#
# ![](https://assets.leetcode.com/uploads/2018/12/14/binarytree.png)
#
#
#
#     Input: root = [3,5,1,6,2,0,8,null,null,7,4], p = 5, q = 4
#     Output: 5
#     Explanation: The LCA of nodes 5 and 4 is 5, since a node can be a descendant of itself according to the LCA definition.
#
#
# **Example 3:**
#
#
#
#     Input: root = [1,2], p = 1, q = 2
#     Output: 1
#
#
#
#
# **Constraints:**
#
#   * The number of nodes in the tree is in the range `[2, 105]`.
#   * `-109 <= Node.val <= 109`
#   * All `Node.val` are **unique**.
#   * `p != q`
#   * `p` and `q` will exist in the tree.
#
#
## @lc code=start
using LeetCode

function lowest_common_ancestor_236(
    root::TreeNode{Int}, p::TreeNode{Int}, q::TreeNode{Int}
)::TreeNode{Int}
    function dfs(node::Union{TreeNode{Int},Nothing})
        (isnothing(node) || node == p || node == q) && return node
        left, right = dfs(node.left), dfs(node.right)
        isnothing(left) && isnothing(right) && return nothing
        !isnothing(left) && !isnothing(right) && return node
        return !isnothing(left) ? left : right
    end
    return dfs(root)
end
## @lc code=end


text:
module LanguageServer
using JSON, REPL, CSTParser, JuliaFormatter, SymbolServer, StaticLint
using CSTParser: EXPR, Tokenize.Tokens, Tokenize.Tokens.kind, headof, parentof, valof, to_codeobject
using StaticLint: refof, scopeof, bindingof
using UUIDs
using Base.Docs, Markdown
import JSONRPC
using JSONRPC: Outbound, @dict_readable

export LanguageServerInstance, runserver

include("URIs2/URIs2.jl")
using .URIs2

JSON.lower(uri::URI) = string(uri)

include("exception_types.jl")
include("protocol/protocol.jl")
include("extensions/extensions.jl")
include("textdocument.jl")
include("document.jl")
include("languageserverinstance.jl")
include("multienv.jl")
include("runserver.jl")
include("staticlint.jl")

include("requests/misc.jl")
include("requests/textdocument.jl")
include("requests/features.jl")
include("requests/hover.jl")
include("requests/completions.jl")
include("requests/workspace.jl")
include("requests/actions.jl")
include("requests/init.jl")
include("requests/signatures.jl")
include("requests/highlight.jl")
include("utilities.jl")

end


text:
include(joinpath("src", "axelrodnet.jl"))

Axelrodnet.run_experiment(
    experiment_name="test_run_grid",
    agentcount=100,
    n_iter=5_000_000,
    nettopology="grid",
    networkprops=Dict("grid_height" => 10),
    stubborncount=0,
    rndseed=0,
    n_replicates=10,
    export_every_n=1_000_000,
    exportdata=true,
    exportmode="default",
    keep_rawnet=false
)

Axelrodnet.run_experiment(
    experiment_name="test_run_grid",
    agentcount=100,
    n_iter=1_000_000,
    nettopology="grid",
    networkprops=Dict("grid_heit" => 10),
    stubborncount=1001,
    rndseed=0,
    n_replicates=10,
    export_every_n=100_000,
    exportdata=true,
    exportmode="aggregated",
    keep_rawnet=true
)

Axelrodnet.run_experiment(
    experiment_name="test_run_watts_strogatz",
    agentcount=100,
    n_iter=1_000_000,
    nettopology="watts_strogatz",
    networkprops=Dict("k" => 101, "beta" => 1.2),
    stubborncount=1001,
    rndseed=0,
    n_replicates=10,
    export_every_n=1_000_000,
    exportdata=true,
    exportmode="default",
    keep_rawnet=true
)

Axelrodnet.run_experiment(
    experiment_name="test_run_barabasi_albert",
    agentcount=100,
    n_iter=1_000_000,
    nettopology="barabasi_albert",
    networkprops=Dict("m0" => 200),
    stubborncount=1001,
    rndseed=0,
    n_replicates=10,
    export_every_n=100_000,
    exportdata=true,
    exportmode="aggregated",
    keep_rawnet=true
)

max_stars_repo_name: JohannesNakayama/axelrodnet

max_issues_repo_name: JohannesNakayama/axelrodnet

max_forks_repo_name: JohannesNakayama/axelrodnet


text:
rng = static_rng()
@test isa(rng, StaticTools.StaticRNG)
@test isa(rng, Xoshiro256✴✴)
@test isa(rand(rng), Float64)
@test 0 <= rand(rng) <= 1


text:
# Creating a test problem
n = 100;
p = 5;
U = randn(n,p);
V = randn(n,p);
K = DiaSymSemiseparable(U,V,ones(n));
x = randn(n);
Kfull = tril(U*V') + triu(V*U',1) + I;

# Testing multiplication
@test isapprox(K*x, Kfull*x, atol = 1e-6)

#  Testing multiplication with the adjoint operator
@test isapprox(K'*x, Kfull'*x, atol = 1e-6)


text:
"""
Description

      The first() method returns the first n number of observations in a TimeArray

Usage

      using MarketData
      first(cl)     # returns the first observation
      first(cl, 10) # returns the first 10 observations

Method Signature(s)

      first(ts::TimeSeries.TimeArray, n::Int=1)

Details

      The first function subsets a TimeSeries.Timearray by the first n observations. When n is not supplied, it
      defaults to n=1.

See Also

      last() which takes the last n observations from a TimeSeries.TimeArray

"""

function first(ts::TimeSeries.TimeArray, n::Int=1)
   ts[1:n]
end


text:
using BDF, Compat.Test

#test BDF read and write
origFilePath = joinpath(dirname(@__FILE__), "Newtest17-256.bdf")
copyFilePath = joinpath(dirname(@__FILE__), "Newtest17-256_copy.bdf")
bdfHeader = readBDFHeader(origFilePath)
dats, evtTab, trigs, statusChan = readBDF(origFilePath)
writeBDF(copyFilePath, dats, trigs, statusChan, bdfHeader["sampRate"][1],
         subjID=bdfHeader["subjID"], recID=bdfHeader["recID"],
         startDate=bdfHeader["startDate"],
         startTime=bdfHeader["startTime"],
         versionDataFormat=bdfHeader["versionDataFormat"],
         chanLabels=bdfHeader["chanLabels"][1:end-1],
         transducer=bdfHeader["transducer"][1:end-1],
         physDim=bdfHeader["physDim"][1:end-1],
         physMin=bdfHeader["physMin"][1:end-1],
         physMax=bdfHeader["physMax"][1:end-1],
         prefilt=bdfHeader["prefilt"][1:end-1],
         reserved=bdfHeader["reserved"])

bdfHeader2 = readBDFHeader(copyFilePath)
dats2, evtTab2, trigs2, statusChan2 = readBDF(copyFilePath)
rm(copyFilePath)
bdfHeader["fileName"] = copyFilePath

@test isequal(bdfHeader2, bdfHeader)
@test isequal(dats2, dats)
@test isequal(evtTab2, evtTab)
@test isequal(trigs2, trigs)
@test isequal(statusChan2, statusChan)


text:
"""
Package: Forecast

    sma(x, n)
    sma(x, n, center)

Smooth a vector of data using a simple moving average.

# Arguments
- `x`: Vector of data.
- `n`: Size of the moving average.
- `center`: if true centers the moving averaged values in the response padding with `missing` values, otherwise the padding takes place at the end.

# Returns
Vector of moving average smoothed values containing `missing` values to preserve the size of the original vector.

# Examples
```julia-repl
julia> sma(1:5,3,true)
5-element Array{Any,1}:
  missing
 2.0
 3.0
 4.0
  missing
```
"""
function sma(x::AbstractVector{T}, n::Integer) where T<:Number

    n == 1 && return x
    N = length(x)
    @assert 1 <= n <= N

    V = Base.promote_op(/, T, typeof(n))
    res = Vector{V}(undef, N-n+1)

    # initial moving average value
    res[1] = ma = sum(x[1:n])/n
    for i in 1:N-n
        res[1+i] = ma = ma + (x[n+i] - x[i]) / n
    end

    res

end

function sma(x::AbstractVector{<:Union{Missing,T}},
             n::Integer,
             center::Bool) where T<:Number

    n == 1 && return x
    N = length(x)
    @assert 1 <= n <= N

    res = sma(collect(skipmissing(x)),n)

    # Missing padding
    ivp = repeat([missing],  findfirst(!ismissing, x)-1 + n÷2)
    fvp = repeat([missing], N-findlast(!ismissing, x)-1 + n-n÷2)

    center ? vcat(ivp,res,fvp) : vcat(res,ivp,fvp)

end


text:
#=
    Implements mpohamiltonian algebra (addition,subtraction,multiplication)
    other actions on mpohamiltonian objects
    the mpohamiltonian transfers
=#

#addition / substraction
function Base.:+(a::MPOHamiltonian{S,T,E},e::AbstractArray{V,1}) where {S,T,E,V}
    length(e) == a.period || throw(ArgumentError("periodicity should match $(a.period) ≠ $(length(e))"))

    nOs = copy(a.Os) # we don't want our addition to change different copies of the original hamiltonian

    for c = 1:a.period
        if nOs[c,1,end] isa E
            nOs[c,1,end] += e[c]
        else
            nOs[c,1,end] += e[c]*isomorphism(storagetype(nOs[c,1,end]),codomain(nOs[c,1,end]),domain(nOs[c,1,end]))
        end
    end

    return MPOHamiltonian{S,T,E}(nOs,a.domspaces,a.pspaces)
end
Base.:-(e::Array{T,1},a::MPOHamiltonian) where T = -1.0*a + e
Base.:+(e::Array{T,1},a::MPOHamiltonian) where T = a + e
Base.:-(a::MPOHamiltonian,e::AbstractArray{T,1}) where T = a + (-e)

function Base.:+(a::MPOHamiltonian{S,T,E},b::MPOHamiltonian{S,T,E}) where {S,T,E}
    a.period == b.period || throw(ArgumentError("periodicity should match $(a.period) ≠ $(b.period)"))
    @assert sanitycheck(a)
    @assert sanitycheck(b)

    #new odim;domspaces,imspaces (pspaces necessarily stays the same)
    nodim = a.odim+b.odim-2
    ndomspaces = PeriodicArray{S,2}(undef,a.period,nodim)

    nOs = PeriodicArray{Union{E,T},3}(fill(zero(E),a.period,nodim,nodim))

    for pos in 1:a.period
        ndomspaces[pos,1:(a.odim-1)] = a.domspaces[pos,1:a.odim-1]
        ndomspaces[pos,a.odim:a.odim+b.odim-2] = b.domspaces[pos,2:b.odim]

        for (i,j) in keys(a,pos)
            #A block
            if(i<a.odim && j<a.odim)
                nOs[pos,i,j] = a.Os[pos,i,j]
            end

            #right side
            if(i<a.odim && j==a.odim)
                nOs[pos,i,nodim] = a.Os[pos,i,j]
            end
        end

        for (i,j) in keys(b,pos)

            #upper Bs
            if(i==1 && j>1)
                if nOs[pos,1,a.odim+j-2] isa T
                    nOs[pos,1,a.odim+j-2] += b[pos,i,j]
                else
                    if b.Os[pos,i,j] isa T

                        #one doesn't exist when dom != codom; hacky workaround
                        #nOs[pos,1,a.odim+j-2] = nOs[pos,1,a.odim+j-2]*one(b[pos,i,j])+b[pos,i,j]
                        if nOs[pos,1,a.odim+j-2] == zero(E)
                            nOs[pos,1,a.odim+j-2] = b[pos,i,j]
                        else
                            t_b = isomorphism(storagetype(b[pos,i,j]),codomain(b[pos,i,j]),domain(b[pos,i,j]));
                            nOs[pos,1,a.odim+j-2] = nOs[pos,1,a.odim+j-2]*t_b+b[pos,i,j]
                        end

                    else
                        nOs[pos,1,a.odim+j-2] += b.Os[pos,i,j]
                    end
                end
            end

            #B block
            if(i>1 && j>1)
                nOs[pos,a.odim+i-2,a.odim+j-2] = b.Os[pos,i,j]
            end
        end
    end

    return MPOHamiltonian{S,T,E}(nOs,ndomspaces,a.pspaces)
end
Base.:-(a::MPOHamiltonian,b::MPOHamiltonian) = a+(-1.0*b)

#multiplication
Base.:*(b::Number,a::MPOHamiltonian)=a*b
function Base.:*(a::MPOHamiltonian{S,T,E},b::Number) where {S,T,E}
    nOs = copy(a.Os)

    for i=1:a.period,j = 1:(a.odim-1)
        nOs[i,j,a.odim]*=b;
    end

    return MPOHamiltonian{S,T,E}(nOs,a.domspaces,a.pspaces)
end

function Base.:*(b::MPOHamiltonian{S,T,E},a::MPOHamiltonian{S,T,E}) where {S,T,E}
    nodim = a.odim*b.odim
    indmap = LinearIndices((a.odim,b.odim))
    nOs = PeriodicArray{Union{E,T},3}(fill(zero(E),a.period,nodim,nodim))

    fusers = PeriodicArray(map(product(1:a.period,1:a.odim,1:b.odim)) do (pos,i,j)
        isomorphism(fuse(a.domspaces[pos,i]*b.domspaces[pos,j]),a.domspaces[pos,i]*b.domspaces[pos,j])
    end)

    ndomspaces = PeriodicArray{S,2}(undef,a.period,nodim)
    for pos = 1:a.period,i in 1:a.odim, j = 1:b.odim
        ndomspaces[pos,indmap[i,j]] = codomain(fusers[pos,i,j])
    end

    for pos = 1:a.period,
        (i,j) in keys(a,pos),
        (k,l) in keys(b,pos)

        if isscal(a,pos,i,j) && isscal(b,pos,k,l)
            nOs[pos,indmap[i,k],indmap[j,l]] = a.Os[pos,i,j]*b.Os[pos,k,l]
        else
            @plansor nOs[pos,indmap[i,k],indmap[j,l]][-1 -2;-3 -4] :=
                fusers[pos,i,k][-1;1 2]*conj(fusers[pos+1,j,l][-4;3 4])*a[pos,i,j][1 5;-3 3]*b[pos,k,l][2 -2;5 4]
        end
    end

    return MPOHamiltonian{S,T,E}(nOs,ndomspaces,a.pspaces)
end

#without the copy, we get side effects when repeating + setindex
Base.repeat(x::MPOHamiltonian{S,T,E},n::Int) where {S,T,E} =
    MPOHamiltonian{S,T,E}(repeat(x.Os,n,1,1),repeat(x.domspaces,n,1),repeat(x.pspaces,n))

function Base.conj(a::MPOHamiltonian)
    b = copy(a.Os)

    for (i,j,k) in keys(a)
        @plansor b[i,j,k][-1 -2;-3 -4]:=conj(a[i,j,k][-1 -3;-2 -4])
    end

    MPOHamiltonian(b)
end

#needed this; perhaps move to tensorkit?
TensorKit.fuse(f::T) where T<: VectorSpace = f


text:
@testset "show" begin
    root = "assets"

    sec = DemoSection(joinpath(root, "section", "default"))
    @test_reference joinpath("references", "section.txt") sec

    page = DemoPage(joinpath(root, "page", "default"))
    sys = Base.Sys.iswindows() ? "windows" : "linux"
    @test_reference joinpath("references", "page_$(sys).txt") page
end


text:
function Dates.format(io, d::Dates.DatePart{'t'}, ep)
    print(io, timescale(ep))
end

abstract type DayOfYearToken end

@inline function Dates.tryparsenext(d::Dates.DatePart{'D'}, str, i, len, locale)
    next = Dates.tryparsenext_base10(str, i, len, 1, 3)
    next === nothing && return nothing
    val, idx = next
    (val >= 1 && val <= 366) || throw(ArgumentError("Day number must be within 1 and 366."))
    return val, idx
end

function Dates.format(io, d::Dates.DatePart{'D'}, ep)
    print(io, dayofyear(ep))
end

abstract type FractionOfSecondToken end

@inline function Dates.tryparsenext(d::Dates.DatePart{'f'}, str, i, len, locale)
    next = Dates.tryparsenext_base10(str, i, len)
    next === nothing && return nothing
    val, idx = next
    val /= 10^(idx - i)
    return val, idx
end

function Dates.format(io, d::Dates.DatePart{'f'}, ep)
    sec = round(Int, fractionofsecond(ep) * 10^d.width)
    print(io, lpad(sec, d.width, '0'))
end

max_stars_repo_name: JuliaAstrodynamics/AstronomicalTime

max_issues_repo_name: JuliaAstrodynamics/AstronomicalTime

max_forks_repo_name: JuliaAstrodynamics/AstronomicalTime


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

include("model_BackendAddressPoolPropertiesFormat.jl")
include("model_FrontendIPConfigurationPropertiesFormat.jl")
include("model_InboundNatPoolPropertiesFormat.jl")
include("model_InboundNatRulePropertiesFormat.jl")
include("model_LoadBalancerListResult.jl")
include("model_LoadBalancerPropertiesFormat.jl")
include("model_LoadBalancingRulePropertiesFormat.jl")
include("model_OutboundNatRulePropertiesFormat.jl")
include("model_ProbePropertiesFormat.jl")
include("model_Resource.jl")
include("model_SubResource.jl")
include("model_BackendAddressPool.jl")
include("model_FrontendIPConfiguration.jl")
include("model_InboundNatPool.jl")
include("model_InboundNatRule.jl")
include("model_LoadBalancer.jl")
include("model_LoadBalancingRule.jl")
include("model_NetworkInterfaceIPConfiguration.jl")
include("model_OutboundNatRule.jl")
include("model_Probe.jl")
include("model_PublicIPAddress.jl")
include("model_Subnet.jl")


text:
using Assignment02
using Test
using Random
using Suppressor

@testset "Assignment02" begin

@testset "Setup" begin
    @test isdefined(Assignment02, :generate_sequence)
    s = generate_sequence(4)
    @test length(s) == 4
    # # Doesn't make sense to do this since rng changes with julia versions
    # Random.seed!(1)
    # @test s == "GGTC"
end

@testset "Question 1" begin
    @test question1 isa DataType
    b = try
        a = question1[]
        push!(a, "test")
    catch e
        println(e)
        a = []
    end
    @test length(b) == 1
end

@testset "Question 2" begin
    @test question2 isa Function
    @test question2("string") isa Bool
    @test question2("string")
    @test question2(1.) isa Bool
    @test !question2(1.)
end

@testset "Question 3" begin
    @test question3 isa Function
    @test question3("AGGC") == 0.75
end

@testset "Question 4" begin
    @test question4 isa Function

    out = @capture_out begin question4("ATTC") end
    (ln1,ln2,ln3,ln4) = split(out, '\n')

    @test ln1 == "Sequence:"
    @test ln2 == "ATTC"
    @test ln3 == "GC Content:"
    @test ln4 == "0.25"

    seq = generate_sequence(20)

    out = @capture_out begin question4(seq) end
    (ln1,ln2,ln3,ln4) = split(out, '\n')

    @test ln1 == "Sequence:"
    @test ln2 == seq
    @test ln3 == "GC Content:"
    @test ln4 == string((count(==('G'), seq)+count(==('C'), seq)) / 20)
end

@testset "Question 5" begin
    @test Assignment02.ce_2_2_1 isa Bool
    @test !Assignment02.ce_2_2_1
    @test Assignment02.ce_2_2_2 isa Bool
    @test Assignment02.ce_2_2_2
    @test Assignment02.ce_2_3_1 isa Float64
    @test isapprox(Assignment02.ce_2_3_1, 4 / 3 * π * 5^3, atol = 1e-2)
    @test Assignment02.ce_2_3_2 isa Float64
    @test isapprox(Assignment02.ce_2_3_2, (24.95 * 0.6) * 60 + 3 + 59 * 0.75, atol = 1e-2)
end

@testset "Question 6" begin
    @test isapprox(bookprice(24.95, 0.4, 60), Assignment02.ce_2_3_2, atol=1e-2)
    @test isapprox(bookprice(24.95, 0.4, 1), 17.97)
    @test isapprox(bookprice(24.95, 0.4, 2), 33.69)
    @test isapprox(bookprice(1, 0, 1), 4.0)
end

end # "Assignment02"