text:
for (F,f) in ((:M,:m), (:S,:s))
    frange = Symbol(f, :range)
    _frange = Symbol(:_, f, :range)
    _fcolon_colon = Symbol(:_, f, :colon_colon)
    fcolon = Symbol(:_, f, :colon)
    _frangestyle = Symbol(:_, f, :rangestyle)
    floatfrange = Symbol(:float, f, :range)
    stepfrangelen_hp = Symbol(:step, f, :rangelen_hp)
    flinspace1 = Symbol(:_, f, :linspace1)
    flinspace = Symbol(:_, f, :linspace)

    LINR = Symbol(:Lin, F, :Range)
    SRL = Symbol(:Step, F, :RangeLen)
    SR = Symbol(:Step, F, :Range)
    UR = Symbol(:Unit, F, :Range)
    @eval begin
        function $(frange)(start; length::Union{Integer,Nothing}=nothing, stop=nothing, step=nothing)
            return $(_frange)(start, step, stop, length)
        end

        function $(frange)(start, stop; length::Union{Integer,Nothing}=nothing, step=nothing)
            return $(_frange)(start, step, stop, length)
        end

        # Range from start to stop: range(a, [step=s,] stop=b), no length
        $(_frange)(start, step,      stop, ::Nothing) = $(fcolon)(start, step, stop)
        $(_frange)(start, ::Nothing, stop, ::Nothing) = $(fcolon)(start, stop)

        $(fcolon)(a::Real, b::Real) = $(fcolon)(promote(a,b)...)

        $(fcolon)(start::T, stop::T) where {T<:Real} = $(UR){T}(start, stop)

        $(fcolon)(start::T, stop::T) where {T} = $(fcolon)(start, oftype(stop-start, 1), stop)

        $(fcolon)(a::T, b::T) where {T<:AbstractFloat} = $(fcolon)(a, T(1), b)

       function $(fcolon)(a::T, b::AbstractFloat, c::T) where {T<:Real}
            return $(fcolon)(promote(a,b,c)...)
        end
        function $(fcolon)(a::T, b::AbstractFloat, c::T) where {T<:AbstractFloat}
            return $(fcolon)(promote(a,b,c)...)
        end
        function $(fcolon)(a::T, b::Real, c::T) where {T<:AbstractFloat}
            return $(fcolon)(promote(a,b,c)...)
        end

        # promote start and stop, leaving step alone
        function $(fcolon)(start::A, step, stop::C) where {A<:Real,C<:Real}
            return $(fcolon)(convert(promote_type(A, C), start),
                             step,
                             convert(promote_type(A, C), stop)
                            )
        end
        function $(fcolon)(start::T, step, stop::T) where {T}
            return $(_fcolon_colon)(start, step, stop)
        end
        function $(fcolon)(start::T, step, stop::T) where {T<:Real}
            return $(_fcolon_colon)(start, step, stop)
        end

        function $(fcolon)(start::T, step::T, stop::T) where T<:Union{Float16,Float32,Float64}
            step == 0 && throw(ArgumentError("range step cannot be zero"))
            # see if the inputs have exact rational approximations (and if so,
            # perform all computations in terms of the rationals)
            step_n, step_d = Base.rat(step)
            if step_d != 0 && T(step_n/step_d) == step
                start_n, start_d = Base.rat(start)
                stop_n, stop_d = Base.rat(stop)
                if start_d != 0 && stop_d != 0 &&
                        T(start_n/start_d) == start && T(stop_n/stop_d) == stop
                    den = lcm(start_d, step_d) # use same denominator for start and step
                    m = maxintfloat(T, Int)
                    if den != 0 && abs(start*den) <= m && abs(step*den) <= m &&  # will round succeed?
                            rem(den, start_d) == 0 && rem(den, step_d) == 0      # check lcm overflow
                        start_n = round(Int, start*den)
                        step_n = round(Int, step*den)
                        len = max(0, div(den*stop_n - stop_d*start_n + step_n*stop_d, step_n*stop_d))
                        # Integer ops could overflow, so check that this makes sense
                        if Base.isbetween(start, start + (len-1)*step, stop + step/2) &&
                                !Base.isbetween(start, start + len*step, stop)
                            # Return a 2x precision range
                            return $(floatfrange)(T, start_n, step_n, len, den)
                        end
                    end
                end
            end
            # Fallback, taking start and step literally
            lf = (stop-start)/step
            if lf < 0
                len = 0
            elseif lf == 0
                len = 1
            else
                len = round(Int, lf) + 1
                stop′ = start + (len-1)*step
                # if we've overshot the end, subtract one:
                len -= (start < stop < stop′) + (start > stop > stop′)
            end
            $(stepfrangelen_hp)(T, start, step, 0, len, 1)
        end
        # AbstractFloat specializations
        function $(_fcolon_colon)(start::T, step, stop::T) where {T}
            T′ = typeof(start+zero(step))
            return $(SR)(convert(T′,start), step, convert(T′,stop))
        end

        function $(fcolon)(start::T, step::T, stop::T) where {T<:AbstractFloat}
            return $(_fcolon_colon)(Base.OrderStyle(T), Base.ArithmeticStyle(T), start, step, stop)
        end
        $(fcolon)(start::T, step::T, stop::T) where {T<:Real} =
            $(_fcolon_colon)(Base.OrderStyle(T), Base.ArithmeticStyle(T), start, step, stop)
        $(_fcolon_colon)(::Base.Ordered, ::Any, start::T, step, stop::T) where {T} = $(SR)(start, step, stop)
        # for T<:Union{Float16,Float32,Float64} see twiceprecision.jl
        $(_fcolon_colon)(::Base.Ordered, ::Base.ArithmeticRounds, start::T, step, stop::T) where {T} =
            $(SRL)(start, step, floor(Int, (stop-start)/step)+1)
        $(_fcolon_colon)(::Any, ::Any, start::T, step, stop::T) where {T} =
            $(SRL)(start, step, floor(Int, (stop-start)/step)+1)
        # range of a given length: range(a, [step=s,] length=l), no stop
        $(_frange)(a::Real,          ::Nothing,         ::Nothing, len::Integer) = $(UR){typeof(a)}(a, oftype(a, a+len-1))
        $(_frange)(a::AbstractFloat, ::Nothing,         ::Nothing, len::Integer) = $(_frange)(a, oftype(a, 1),   nothing, len)
        $(_frange)(a::AbstractFloat, st::AbstractFloat, ::Nothing, len::Integer) = $(_frange)(promote(a, st)..., nothing, len)
        $(_frange)(a::Real,          st::AbstractFloat, ::Nothing, len::Integer) = $(_frange)(float(a), st,      nothing, len)
        $(_frange)(a::AbstractFloat, st::Real,          ::Nothing, len::Integer) = $(_frange)(a, float(st),      nothing, len)
        $(_frange)(a,                ::Nothing,         ::Nothing, len::Integer) = $(_frange)(a, oftype(a-a, 1), nothing, len)

        function $(_frange)(a::T, step::T, ::Nothing, len::Integer) where {T <: AbstractFloat}
            return $(_frangestyle)(Base.OrderStyle(T), Base.ArithmeticStyle(T), a, step, len)
        end
        function $(_frange)(a::T, step, ::Nothing, len::Integer) where {T}
            return $(_frangestyle)(Base.OrderStyle(T), Base.ArithmeticStyle(T), a, step, len)
        end
        function $(_frangestyle)(::Base.Ordered, ::Base.ArithmeticWraps, a::T, step::S, len::Integer) where {T,S}
            return $(SR){T,S}(a, step, convert(T, a+step*(len-1)))
        end
        function $(_frangestyle)(::Any, ::Any, a::T, step::S, len::Integer) where {T,S}
            return $(SRL){typeof(a+0*step),T,S}(a, step, len)
        end

        function $(_frange)(start::T, ::Nothing, stop::T, len::Integer) where {T<:IEEEFloat}
            len < 2 && return $(flinspace1)(T, start, stop, len)
            start == stop && $(stepfrangelen_hp)(T, start, zero(T), 0, len, 1)
            # Attempt to find exact rational approximations
            start_n, start_d = rat(start)
            stop_n, stop_d = rat(stop)
            if start_d != 0 && stop_d != 0
                den = lcm(start_d, stop_d)
                m = maxintfloat(T, Int)
                if den != 0 && abs(den*start) <= m && abs(den*stop) <= m
                    start_n = round(Int, den*start)
                    stop_n = round(Int, den*stop)
                    if T(start_n/den) == start && T(stop_n/den) == stop
                        return $(flinspace)(T, start_n, stop_n, len, den)
                    end
                end
            end
            return $(flinspace)(start, stop, len)
        end

        $(_frange)(start,     step,      ::Nothing, ::Nothing) = # range(a, step=s)
            throw(ArgumentError("At least one of `length` or `stop` must be specified"))
        $(_frange)(start,     ::Nothing, ::Nothing, ::Nothing) = # range(a)
            throw(ArgumentError("At least one of `length` or `stop` must be specified"))
        $(_frange)(::Nothing, ::Nothing, ::Nothing, ::Nothing) = # range(nothing)
            throw(ArgumentError("At least one of `length` or `stop` must be specified"))
        $(_frange)(start::Real, step::Real, stop::Real, length::Integer) = # range(a, step=s, stop=b, length=l)
            throw(ArgumentError("Too many arguments specified; try passing only one of `stop` or `length`"))
        $(_frange)(::Nothing, ::Nothing, ::Nothing, ::Integer) = # range(nothing, length=l)
            throw(ArgumentError("Can't start a range at `nothing`"))

        function $(_frange)(start::T, ::Nothing, stop::S, len::Integer) where {T,S}
            a, b = promote(start, stop)
            return $(_frange)(a, nothing, b, len)
        end
        $(_frange)(start::T, ::Nothing, stop::T, len::Integer) where {T<:Real} = $(LINR){T}(start, stop, len)
        $(_frange)(start::T, ::Nothing, stop::T, len::Integer) where {T} = $(LINR){T}(start, stop, len)
        $(_frange)(start::T, ::Nothing, stop::T, len::Integer) where {T<:Integer} =
            $(flinspace)(float(T), start, stop, len)
        ## for Float16, Float32, and Float64 we hit twiceprecision.jl to lift to higher precision StepRangeLen
        # for all other types we fall back to a plain old LinRange
        $(flinspace)(::Type{T}, start::Integer, stop::Integer, len::Integer) where T = $(LINR){T}(start, stop, len)

        function $(flinspace)(::Type{T}, start::Integer, stop::Integer, len::Integer) where {T<:IEEEFloat}
            return $(flinspace)(T, start, stop, len, 1)
        end
        function $(flinspace)(
            ::Type{T},
            start_n::Integer,
            stop_n::Integer,
            len::Integer,
            den::Integer
           ) where {T<:IEEEFloat}
            len < 2 && return $(flinspace1)(T, start_n/den, stop_n/den, len)
            start_n == stop_n && return $(stepfrangelen_hp)(T, (start_n, den), (zero(start_n), den), 0, len, 1)
            tmin = -start_n/(Float64(stop_n) - Float64(start_n))
            imin = round(Int, tmin*(len-1)+1)
            imin = clamp(imin, 1, Int(len))
            ref_num = Int128(len-imin) * start_n + Int128(imin-1) * stop_n
            ref_denom = Int128(len-1) * den
            ref = (ref_num, ref_denom)
            step_full = (Int128(stop_n) - Int128(start_n), ref_denom)
            return $(stepfrangelen_hp)(T, ref, step_full,  nbitslen(T, len, imin), Int(len), imin)
        end

        function $(flinspace1)(::Type{T}, start, stop, len::Integer) where {T<:IEEEFloat}
            len >= 0 || throw(ArgumentError("range($start, stop=$stop, length=$len): negative length"))
            if len <= 1
                len == 1 && (start == stop || throw(ArgumentError("range($start, stop=$stop, length=$len): endpoints differ")))
                # Ensure that first(r)==start and last(r)==stop even for len==0
                # The output type must be consistent with steprangelen_hp
                if T<:Union{Float32,Float16}
                    return $(SRL){T}(Float64(start), Float64(start) - Float64(stop), len, 1)
                else
                    return $(SRL)(TwicePrecision(start, zero(T)), TwicePrecision(start, -stop), len, 1)
                end
            end
            throw(ArgumentError("should only be called for len < 2, got $len"))
        end
        function $(_frange)(a::T, st::T, ::Nothing, len::Integer) where T<:Union{Float16,Float32,Float64}
            start_n, start_d = Base.rat(a)
            step_n, step_d = Base.rat(st)
            if start_d != 0 && step_d != 0 &&
                    T(start_n/start_d) == a && T(step_n/step_d) == st
                den = lcm(start_d, step_d)
                m = Base.maxintfloat(T, Int)
                if abs(den*a) <= m && abs(den*st) <= m &&
                        rem(den, start_d) == 0 && rem(den, step_d) == 0
                    start_n = round(Int, den*a)
                    step_n = round(Int, den*st)
                    return $(floatfrange)(T, start_n, step_n, len, den)
                end
            end
            return $(stepfrangelen_hp)(T, a, st, 0, len, 1)
        end

        function $(stepfrangelen_hp)(
            ::Type{Float64},
            ref::Tuple{Integer,Integer},
            step::Tuple{Integer,Integer},
            nb::Integer,
            len::Integer,
            offset::Integer
           )
            $(SRL)(TwicePrecision{Float64}(ref), TwicePrecision{Float64}(step, nb), Int(len), offset)
        end

        function $(stepfrangelen_hp)(
            ::Type{T},
            ref::Tuple{Integer,Integer},
            step::Tuple{Integer,Integer},
            nb::Integer,
            len::Integer,
            offset::Integer
           ) where {T<:IEEEFloat}
            return $(SRL){T}(ref[1]/ref[2], step[1]/step[2], Int(len), offset)
        end

        function $(stepfrangelen_hp)(
            ::Type{Float64},
            ref::Base.F_or_FF,
            step::Base.F_or_FF,
            nb::Integer,
            len::Integer,
            offset::Integer
           )
            return $(SRL)(Base._TP(ref), twiceprecision(Base._TP(step), nb), Int(len), offset)
        end

        function $(stepfrangelen_hp)(
            ::Type{T}, ref::Base.F_or_FF,
            step::Base.F_or_FF, nb::Integer,
            len::Integer,
            offset::Integer
           ) where {T<:IEEEFloat}
            return $(SRL){T}(Base.asF64(ref), Base.asF64(step), Int(len), offset)
        end
        function $(floatfrange)(::Type{T}, start_n::Integer, step_n::Integer, len::Integer, den::Integer) where T
            if len < 2 || step_n == 0
                return $(stepfrangelen_hp)(T, (start_n, den), (step_n, den), 0, Int(len), 1)
            end
            # index of smallest-magnitude value
            imin = clamp(round(Int, -start_n/step_n+1), 1, Int(len))
            # Compute smallest-magnitude element to 2x precision
            ref_n = start_n+(imin-1)*step_n  # this shouldn't overflow, so don't check
            nb = nbitslen(T, len, imin)
            return $(stepfrangelen_hp)(T, (ref_n, den), (step_n, den), nb, Int(len), imin)
        end

        function $(floatfrange)(a::AbstractFloat, st::AbstractFloat, len::Real, divisor::AbstractFloat)
            T = promote_type(typeof(a), typeof(st), typeof(divisor))
            m = maxintfloat(T, Int)
            if abs(a) <= m && abs(st) <= m && abs(divisor) <= m
                ia, ist, idivisor = round(Int, a), round(Int, st), round(Int, divisor)
                if ia == a && ist == st && idivisor == divisor
                    # We can return the high-precision range
                    return $(floatfrange)(T, ia, ist, Int(len), idivisor)
                end
            end
            # Fallback (misses the opportunity to set offset different from 1,
            # but otherwise this is still high-precision)
            return $(stepfrangelen_hp)(T, (a,divisor), (st,divisor), Base.nbitslen(T, len, 1), Int(len), 1)
        end
        function $(flinspace)(start::T, stop::T, len::Integer) where {T<:IEEEFloat}
            (isfinite(start) && isfinite(stop)) || throw(ArgumentError("start and stop must be finite, got $start and $stop"))
            # Find the index that returns the smallest-magnitude element
            Δ, Δfac = stop-start, 1
            if !isfinite(Δ)   # handle overflow for large endpoints
                Δ, Δfac = stop/len - start/len, Int(len)
            end
            tmin = -(start/Δ)/Δfac            # t such that (1-t)*start + t*stop == 0
            imin = round(Int, tmin*(len-1)+1) # index approximately corresponding to t
            if 1 < imin < len
                # The smallest-magnitude element is in the interior
                t = (imin-1)/(len-1)
                ref = T((1-t)*start + t*stop)
                step = imin-1 < len-imin ? (ref-start)/(imin-1) : (stop-ref)/(len-imin)
            elseif imin <= 1
                imin = 1
                ref = start
                step = (Δ/(len-1))*Δfac
            else
                imin = Int(len)
                ref = stop
                step = (Δ/(len-1))*Δfac
            end
            if len == 2 && !isfinite(step)
                # For very large endpoints where step overflows, exploit the
                # split-representation to handle the overflow
                return $(stepfrangelen_hp)(T, start, (-start, stop), 0, 2, 1)
            end
            # 2x calculations to get high precision endpoint matching while also
            # preventing overflow in ref_hi+(i-offset)*step_hi
            m, k = prevfloat(floatmax(T)), max(imin-1, len-imin)
            step_hi_pre = clamp(step, max(-(m+ref)/k, (-m+ref)/k), min((m-ref)/k, (m+ref)/k))
            nb = nbitslen(T, len, imin)
            step_hi = truncbits(step_hi_pre, nb)
            x1_hi, x1_lo = add12((1-imin)*step_hi, ref)
            x2_hi, x2_lo = add12((len-imin)*step_hi, ref)
            a, b = (start - x1_hi) - x1_lo, (stop - x2_hi) - x2_lo
            step_lo = (b - a)/(len - 1)
            ref_lo = a - (1 - imin)*step_lo
            return $(stepfrangelen_hp)(T, (ref, ref_lo), (step_hi, step_lo), 0, Int(len), imin)
        end
    end
end


text:
# Copyright 2020 Google LLC

# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at

#     https://www.apache.org/licenses/LICENSE-2.0

# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.

input = readlines("day01input.txt")
nums = map(s->parse(Int, s), input)

function part1(nums)
    for i=1:length(nums)-1
        for j=(i+1):length(nums)
            if nums[i] + nums[j] == 2020
                return nums[i] * nums[j]
            end
        end
    end
end

function part2(nums)
    for i=1:length(nums)-2
        for j=(i+1):length(nums)-1
            for k=(j+1):length(nums)
                if nums[i] + nums[j] + nums[k] == 2020
                    return nums[i] * nums[j] * nums[k]
                end
            end
        end
    end
end

part1(nums)
part2(nums)

ext: jl


text:
function test_fibonacci_nthfibonacci()
    @test nthfibonacci(1) == 1
    @test nthfibonacci(2) == 1
    @test nthfibonacci(10) == 55

    @test nthfibonacci(91, Int64) == 4660046610375530309
    @test_throws OverflowError nthfibonacci(92, Int64)

    @test nthfibonacci(183, Int128) == 78569350599398894027251472817058687522
    @test_throws OverflowError nthfibonacci(184, Int128)

    @test nthfibonacci(200, BigInt) == 280571172992510140037611932413038677189525
end

function test_fibonacci_nfibonacci()
    @test nfibonacci(10) == [1, 1, 2, 3, 5, 8, 13, 21, 34, 55]
    @test eltype(nfibonacci(10)) == typeof(1)
end

function test_fibonacci_somefibonacci()
    @test collect(somefibonacci(34)) == [1, 1, 2, 3, 5, 8, 13, 21, 34]
    @test collect(somefibonacci(39, 13, -5)) == [13, -5, 8, 3, 11, 14, 25, 39]
    @test eltype(somefibonacci(34)) == typeof(1)
    @test eltype(somefibonacci(39, 13, -5)) == typeof(13)
end

function test_fibonacci_exactfibonacci()
    @test collect(exactfibonacci(5)) == [1, 1, 2, 3, 5]
    @test collect(exactfibonacci(5, 13, -5)) == [13, -5, 8, 3, 11]
    @test eltype(exactfibonacci(5)) == typeof(1)
    @test eltype(exactfibonacci(5, 13, -5)) == typeof(13)
end

function test_fibonacci_all()
    print(rpad("Math.NumberTheory.Fibonacci...", 50, ' '))

    test_fibonacci_nthfibonacci()
    test_fibonacci_nfibonacci()
    test_fibonacci_somefibonacci()
    test_fibonacci_exactfibonacci()

    println("PASS")
end


text:
module LoggingFormats

import Logging, JSON3, StructTypes

###############
## Truncated ##
###############

shorten_str(str, max_len) = shorten_str(string(str), max_len)
function shorten_str(str::String, max_len)
    if textwidth(str) <= max_len
        return SubString(str)
    end
    len = textwidth('…')
    ind = 1
    for i in eachindex(str)
        c = @inbounds str[i]
        len += textwidth(c)
        len > max_len && break
        ind = i
    end
    return SubString(str, 1, ind) * '…'
end

struct Truncated <: Function
    max_var_len::Int
    Truncated(max_var_len) = max_var_len <= 0 ? error("max_var_len must be positive") : new(max_var_len)
end
Truncated() = Truncated(5_000)

# copied from https://github.com/JuliaLang/julia/blob/v1.5.4/stdlib/Logging/src/ConsoleLogger.jl and modified
function (tr::Truncated)(io, args)
    levelstr = args.level == Logging.Warn ? "Warning" : string(args.level)
    msglines = split(chomp(shorten_str(args.message, tr.max_var_len)), '\n')
    println(io, "┌ ", levelstr, ": ", msglines[1])
    for i in 2:length(msglines)
        str_line = sprint(print, "│ ", msglines[i])
        println(io, shorten_str(str_line, tr.max_var_len))
    end
    for (key, val) in args.kwargs
        str_line = sprint(print, "│   ", key, " = ", val)
        println(io, shorten_str(str_line, tr.max_var_len))
    end
    println(io, "└ @ ", something(args._module, "nothing"), " ",
            something(args.file, "nothing"), ":", something(args.line, "nothing"))
    nothing
end

##########
## JSON ##
##########

lvlstr(lvl::Logging.LogLevel) = lvl >= Logging.Error ? "error" :
                                lvl >= Logging.Warn  ? "warn"  :
                                lvl >= Logging.Info  ? "info"  :
                                                       "debug"

struct JSONLogMessage
    level::String
    msg::String
    _module::Union{String,Nothing}
    file::Union{String,Nothing}
    line::Union{Int,Nothing}
    group::Union{String,Nothing}
    id::Union{String,Nothing}
    kwargs::Dict{String,String}
end
function JSONLogMessage(args)
    JSONLogMessage(
        lvlstr(args.level),
        args.message isa AbstractString ? args.message : string(args.message),
        args._module === nothing ? nothing : string(args._module),
        args.file,
        args.line,
        args.group === nothing ? nothing : string(args.group),
        args.id === nothing ? nothing : string(args.id),
        Dict{String,String}(string(k) => string(v) for (k, v) in args.kwargs)
    )
end
StructTypes.StructType(::Type{JSONLogMessage}) = StructTypes.OrderedStruct()
StructTypes.names(::Type{JSONLogMessage}) = ((:_module, :module), )

struct JSON <: Function
end

function (::JSON)(io, args)
    logmsg = JSONLogMessage(args)
    JSON3.write(io, logmsg)
    println(io)
    return nothing
end

############
## logfmt ##
############
# See  https://brandur.org/logfmt

struct LogFmt <: Function
end
function (::LogFmt)(io, args)
    print(io, "level=", lvlstr(args.level),
              " msg=\"",
    )
    escape_string(io, args.message isa AbstractString ? args.message : string(args.message), '"')
    print(io, "\"",
              " module=", something(args._module, "nothing"),
              " file=\"",
    )
    escape_string(io, args.file isa AbstractString ? args.file : string(something(args.file, "nothing")), '"')
    print(io, "\"",
              " line=", something(args.line, "nothing"),
              " group=\"",
    )
    escape_string(io, args.group isa AbstractString ? args.group : string(something(args.group, "nothing")), '"')
    print(io, "\"",
              " id=", something(args.id, "nothing"),
    )
    for (k, v) in args.kwargs
        print(io, " ", k, "=\"")
        escape_string(io, sprint(print, something(v, "nothing")), '"')
        print(io, "\"")
    end
    println(io)
    return nothing
end

end # module


text:
using CSV, DataFrames, ScikitLearn, PyPlot
pathtodata = joinpath("julia-scripts", "model-zoo", "covid_cleaned.csv")
data = DataFrame(CSV.File(pathtodata))

X = convert(Array, data[!, Not(:covid_res)])
y = convert(Array, data[!, :covid_res])

@sk_import model_selection:train_test_split
X_train, X_test, y_train, y_test = train_test_split(X, y; test_size=0.5, random_state=42)

@sk_import naive_bayes:GaussianNB
gnb = GaussianNB()
fit!(gnb, X_train, y_train)

y_predict = predict(gnb, X_train)

@sk_import metrics:classification_report
print(classification_report(y_train, y_predict))

@sk_import metrics:plot_confusion_matrix
plot_confusion_matrix(gnb, X_train, y_train)
PyPlot.gcf()

max_stars_repo_name: coinslab/ComputationalCognitiveModeling

max_issues_repo_name: coinslab/ComputationalCognitiveModeling

max_forks_repo_name: coinslab/ComputationalCognitiveModeling


text:
using Combinatorics

function solve(arr)
	pairs = collect(combinations(arr, 2));
	idx = findfirst(x -> x[1] + x[2] == 2020, pairs);
	pairs[idx][1] * pairs[idx][2]
end

function solve_for_2_better(input)
	arr = sort(input)
	lower_i = 1
	upper_i = length(arr)
	sum = arr[lower_i] + arr[upper_i]
	while sum != 2020
		if sum < 2020
			lower_i += 1
		else
			upper_i -= 1
		end
		sum = arr[lower_i] + arr[upper_i]
	end
	arr[lower_i] * arr[upper_i]
end

function solve_for_2_also_better(input)
	arr = zeros(2020)
	for x in input
		arr[x] = true
		if arr[2020 - x] == true
			return x * (2020 - x)
		end
	end
end

function solve_for_n(arr, n)
	pairs = collect(combinations(arr, n));
	idx = findfirst(x -> sum(x) == 2020, pairs);
	reduce(*, pairs[idx])
end

test_input = [1721, 979, 366, 299, 675, 1456];
@assert solve(test_input) == 514579 ["test case 1 failed"];
@assert solve_for_2_better(test_input) == 514579 ["test case 1a failed"];
@assert solve_for_2_also_better(test_input) == 514579 ["test case 1b failed"];
@assert solve_for_n(test_input, 3) == 241861950 ["test case 2 failed"];

input = map(x -> parse(Int64, x), readlines("input/day1.txt"));
println("First star: " * string(@time solve(input)))
println("First star: " * string(@time solve_for_2_better(input)))
println("First star: " * string(@time solve_for_2_also_better(input)))
println("Second star: " * string(@time solve_for_n(input, 3)))


text:
@inline function _remove_duplicate_spaces(str)::String
    a = replace(str, isspace => " ")
    while occursin("  ", a)
        a = replace(a, "  " => " ")
    end
    return a
end

@inline function _add_space_before_uppercase_characters(c::Char)::String
    if isuppercase(c)
        return string(" ", lowercase(c))
    else
        return string(c)
    end
end

@inline function _add_space_before_uppercase_characters(str::AbstractString)::String
    original_chars = collect(str)
    new_elements = _add_space_before_uppercase_characters.(original_chars)
    result = join(new_elements)
    return result
end

@inline function spacecase(str)::String
    a = strip(str)
    b = replace(str, "-" => " ")
    c = replace(b,   "_" => " ")
    d = replace(c, isspace => " ")
    e = _add_space_before_uppercase_characters(d)
    f = lowercase(e)
    g = _remove_duplicate_spaces(f)
    h = strip(g)
    return h
end

"""
Convert the input string to snake_case.
"""
@inline function snakecase(str)::String
    a = spacecase(str)
    b = replace(a, isspace => "_")
    return b
end

"""
Convert the input string to PascalCase.
"""
@inline function pascalcase(str)::String
    a = spacecase(str)
    b = titlecase(a; strict = true, wordsep = isspace)
    c = replace(b, isspace => "")
    return c
end


text:
using Documenter
using Graphics

makedocs(
    sitename = "Graphics",
    format = Documenter.HTML(prettyurls = get(ENV, "CI", nothing) == "true"),
    modules = [Graphics],
    pages = ["index.md", "reference.md"]
)

deploydocs(
    repo = "github.com/JuliaGraphics/Graphics.jl.git"
)


text:
# filled contours (contourf)

using CairoMakie

xs = LinRange(0, 10, 100)
ys = LinRange(0, 10, 100)
zs = [cos(x) * sin(y) for x in xs, y in ys]

f = Figure()
Axis(f[1, 1])

co = CairoMakie.contourf!(xs, ys, zs, levels = 10, colormap = :flag_ar)

Colorbar(f[1, 2], co)

max_stars_repo_name: Ellon-M/visualizations

max_issues_repo_name: Ellon-M/visualizations

max_forks_repo_name: Ellon-M/visualizations


text:
# This file is a part of SolidStateDetectors.jl, licensed under the MIT License (MIT).

include("Planes.jl")

include("plot_recipes.jl")


text:
"""
    getconstant(s::Union{MOI.EqualTo, MOI.GreaterThan, MOI.LessThan})

Returns the constant of the set.
"""
getconstant(s::MOI.EqualTo) = s.value
getconstant(s::MOI.GreaterThan) = s.lower
getconstant(s::MOI.LessThan) = s.upper

const DimensionUpdatableSets = Union{MOI.Reals,
                                     MOI.Zeros,
                                     MOI.Nonnegatives,
                                     MOI.Nonpositives}
"""
    updatedimension(s::AbstractVectorSet, newdim)

Returns a set with the dimension modified to `newdim`.
"""
function updatedimension(::S, newdim) where S<:DimensionUpdatableSets
    S(newdim)
end

"""
    shift_constant(set::MOI.AbstractScalarSet,
                   offset)

Returns a new scalar set `new_set` such that `func`-in-`set` is equivalent to
`func + offset`-in-`new_set`.

## Examples

The call `shift_constant(MOI.Interval(-2, 3), 1)` is equal to
`MOI.Interval(-1, 4)`.
"""
function shift_constant(set::Union{MOI.LessThan{T},
                                   MOI.GreaterThan{T},
                                   MOI.EqualTo{T}},
                        offset::T) where T
    return typeof(set)(getconstant(set) + offset)
end
function shift_constant(set::MOI.Interval, offset)
    return MOI.Interval(set.lower + offset, set.upper + offset)
end


text:
function steadystate!(m::RealBond;
                      βlo::Float64 = 0.4/m[:R],
                      βhi::Float64 = 1 ./m[:R],
                      excess::Float64 = 5000.,
                      tol::Float64 = 1e-5,
                      maxit::Int64 = 100)

    # Load in parameters and grids
    R     = m[:R]
    γ     = m[:γ]
    ν     = m[:ν]
    abar  = m[:abar]
    sigs  = m[:σ_s]

    elo   = get_setting(m, :elo)
    ehi   = get_setting(m, :ehi)

    xgrid = m.grids[:xgrid].points
    xwts  = m.grids[:xgrid].weights
    sgrid = m.grids[:sgrid].points
    swts  = m.grids[:sgrid].weights
    ggrid = m.grids[:ggrid]

    xgrid_total   = m.grids[:xgrid_total]
    sgrid_total   = m.grids[:sgrid_total]
    weights_total = m.grids[:weights_total]

    nx      = get_setting(m, :nx)
    ns      = get_setting(m, :ns)
    n       = get_setting(m, :n)

    # Flag allows one not to recompute steady state
    if get_setting(m, :use_last_βstar)
       βlo = βhi = m[:βstar]
    end

    # Initialization
    count   = 1
    c = zeros(n)
    η = zeros(n)
    ap = zeros(n)
    KF = zeros(n,n)
    μ = zeros(n)
    β = (βlo + βhi)/2.0

    # Define the mollifier
    qfunction(x) = mollifier_realbond(x, ehi, elo)

    # compute constrained consumption and labor supply once and for all
    χss = zeros(nx*ns)
    # ηss = zeros(nx*ns)
    aborrow  = abar/R

    laborzero(s::AbstractFloat, x::AbstractFloat, c::AbstractFloat,
              ν::AbstractFloat, γ::AbstractFloat, aborrow::AbstractFloat) = s^(1+ν) - ((aborrow+c-x)^ν)*c^γ
    chifun(s::AbstractFloat, x::AbstractFloat, ν::AbstractFloat,
           γ::AbstractFloat, aborrow::AbstractFloat) = fzero(c ->laborzero(s, x, c, ν, γ, aborrow), x-aborrow+1e-8, 15.0)
    # etafun(s::AbstractFloat,x::AbstractFloat,ν::AbstractFloat,γ::AbstractFloat,aborrow::AbstractFloat) = chifun(s,x,ν,γ,aborrow)^(-(γ/ν))*s^(1.0/ν)

    for iss in 1:ns
        for ix in 1:nx
            χss[ix + nx*(iss-1)] = chifun(sgrid[iss],xgrid[ix],ν.value,γ.value,aborrow)
            # ηss[ix + nx*(iss-1)] = etafun(sgrid[iss],xgrid[ix],ν,γ,aborrow)
        end
    end

    # guess for the W function W(w) = β R E u'(c_{t+1})
    βguess = 0.8

    # MDχ: This l_in doesn't seem to be used, since it gets immediately overwritten by
    # l_in from policy, unless the while loop condition is never entered
    l_in_const = βguess*R*sum((qfunction.(abar .- xgrid_total).*χss.^(-γ)).*(kron((swts.*ggrid),xwts)))
    l_in = l_in_const*ones(nx*ns)

    while abs(excess) > tol && count < maxit # clearing markets
#        @warn "You're recomputing steady state! BAD!"
        l_in_guess = ones(n)
        β = (βlo + βhi)/2.0

        (c, η, ap, l_in, KF) = policy_realbond(β, R.value, γ.value, ν.value, sigs.value,
                                               ehi, elo, xgrid_total, sgrid_total,
                                               xwts, swts, l_in_guess, ggrid, χss)

        # sspolicy returns the consumption and hours decision rules, a prime
        # decision rule, l_in = updated β R u'(c_{t+1}),
        # KF is the Kolmogorov foward operator
        # and is the map which moves you from cash on hand distribution
        # today to cash on had dist tomorrow

        LPMKF = xwts[1]*swts[1]*KF

        # find eigenvalue closest to 1
        (Dee,Vee) = eigen(LPMKF)
        if abs(Dee[1]-1)>2e-1 # that's the tolerance we are allowing
            @warn "your eigenvalue is ", Dee[1], " which is too far from 1, something is wrong"
        end

        μ = real(Vee[:,1]) # Pick the eigen vector associated with the largest
                           # eigenvalue and move it back to values
                           # mdχ: the μ chosen is the eigenvector associated to the largest
                           # eigenvalue of the KF because that means, we have induced the
                           # shift in μ by the largest amount possible for a given KF.
                           # Hence, to find the steady-state μ, we want to pick the
                           # potential μ that changes by the most (i.e. the largest
                           # eigenvector) so we have an upper bound on the degree of change
                           # of μ induced by the KF.

        μ = μ/(weights_total'*μ) # Scale of eigenvectors not determinate: rescale to integrate to exactly 1

        excess = weights_total'*(μ.*ap)  # compute excess supply of savings, which is a fn of w

        # bisection
        if excess > 0
            βhi = β
        elseif excess < 0
            βlo = β
        end
        count += 1
    end

    m[:lstar] = l_in
    m[:cstar] = c
    m[:ηstar] = η
    m[:μstar] = μ
    m[:βstar] = β
    m[:χstar] = χss

    nothing
end

function policy_realbond(β::AbstractFloat, R::AbstractFloat,
                         γ::AbstractFloat, ν::AbstractFloat,
                         sigs::AbstractFloat, ehi::AbstractFloat, elo::AbstractFloat,
                         xgrid_total::Vector{Float64}, sgrid_total::Vector{Float64},
                         xwts::Vector{Float64}, swts::Vector{Float64},
                         l_in::Vector{Float64}, ggrid::Vector{Float64},
                         χss::Vector{Float64};
                         damp::Float64 = 0.4, dist::Float64 = 1.,
                         tol::Float64 = 1e-4, maxit::Int64 = 500)
    # outputs: [c, η, ap, l_out, tr]

    nx      = length(xwts)
    ns      = length(swts)
    n       = ns*nx
    c       = zeros(n)      # consumption
    η       = zeros(n)      # hours
    ap      = zeros(n)      # savings
    counter = 1
    l_in    = ones(n)
    l_out   = copy(l_in)    # initialize l_out
    qfunction(x) = mollifier_realbond(x,ehi,elo)

    while dist > tol && counter < maxit
        # compute c(w) given guess for l_in = β*R*E[u'(c_{t+1})]

        for i in 1:n
            c[i]  = min(l_in[i]^(-1.0/γ), χss[i])
        end

        η  = (sgrid_total.^(1.0/ν)).*(c.^(-γ/ν))
        ap = R*(xgrid_total+sgrid_total.*η - c)

        l_out = parameterized_expectations_realbond(nx, ns, qfunction, xgrid_total, ggrid, xwts, swts, ap, c, γ, β, R)

        dist = norm(l_out - l_in, Inf)

        # if counter%1==0
        #     println(dist)
        # end
        l_in = damp*l_out + (1.0-damp)*l_in
        counter += 1
    end

    if counter == maxit
        @warn "Euler iteration did not converge"
    end

    # c  = min.(l_in.^(-1.0/γ),χss)
    # η  = (sgrid_total.^(1.0/ν)).*(c.^(-γ/ν))
    # ap = R*(xgrid_total+sgrid_total.*η - c)

    tr = kolmogorov_fwd_policy_realbond(nx, ns, qfunction, xgrid_total, ggrid, ap)

    return c, η, ap, l_out, tr
end

function parameterized_expectations_realbond(nx::Int, ns::Int,
                                             qfunction::Function,
                                             xgrid_total::Vector{Float64}, ggrid::Vector{Float64},
                                             xwts::Vector{Float64}, swts::Vector{Float64},
                                             ap::Vector{Float64}, c::Vector{Float64},
                                             γ::Float64, β::Float64, R::Float64)
    l_out   = ones(ns*nx)    # initialize l_out
    for ix in 1:nx
        for iss in 1:ns
            sumn = 0.0
            for ixp in 1:nx
                for isp in 1:ns
                    # This is the parameterized expectations algorithm approximation of the
                    # conditional expectation within the Euler equation
                   sumn += qfunction(ap[ix+nx*(iss-1)] - xgrid_total[nx*(isp-1)+ixp])*ggrid[isp]*xwts[ixp]*swts[isp]*(c[nx*(isp-1)+ixp]^(-γ))
                end
            end
            l_out[ix+nx*(iss-1)] = β*R*sumn
        end
    end

    return l_out
end

function kolmogorov_fwd_policy_realbond(nx::Int, ns::Int,
                                        qfunction::Function,
                                        xgrid_total::Vector{Float64}, ggrid::Vector{Float64},
                                        ap::Vector{Float64})
    tr = zeros(ns*nx,ns*nx)    # transition matrix (KFE) mapping todays dist of w to w' (i.e. μ to μ′)
    for ix in 1:nx
        for ixp in 1:nx
            for iss in 1:ns
                for isp in 1:ns
                    tr[ixp+nx*(isp-1),ix+nx*(iss-1)] += qfunction(ap[ix+nx*(iss-1)] - xgrid_total[nx*(isp-1)+ixp])*ggrid[isp]
                end
            end
        end
    end
    return tr
end

function mollifier_realbond(z::AbstractFloat,ehi::AbstractFloat,elo::AbstractFloat)
    # mollifier function
    In = 0.443993816237631
    if z<ehi && z>elo
        temp = -1.0 + 2.0*(z-elo)/(ehi-elo)
        out = (2.0/(ehi-elo))*exp(-1.0/(1.0-temp*temp))/In
    else
        out = 0.0
    end
    return out
end

function dmollifier_realbond(x::AbstractFloat, ehi::AbstractFloat, elo::AbstractFloat)
    In = 0.443993816237631
    if x<ehi && x>elo
        temp = (-1.0 + 2.0*(x-elo)/(ehi-elo))
        out  = -(2*temp./((1 - temp.^2).^2)).*(2/(ehi-elo)).*mollifier_realbond(x, ehi, elo)
    else
        out = 0.0
    end
    return out
end


text:
module TestDoctest

import BangBang
using Documenter: doctest
using Test

@testset "doctest" begin
    doctest(BangBang; manual = true)
end

end  # module


text:
using Test, Plots, LinearAlgebra, Revise, FourierGPE

# Test incompressible without loading VortexDistributions
abstract type Vortex end

struct PointVortex <: Vortex
    xv::Float64
    yv::Float64
    qv::Int64
end

Λ = 0.8249
r(x,y) = sqrt(x^2+y^2)
vcore(r) = (Λ*r)^2/(1+(Λ*r)^2)
vdensity(x,y) = vcore(r(x,y))
vphase(x,y) = atan(y,x)
vortexgrid(x,y,x0,y0,q0) = vdensity(x - x0, y - y0)*exp(im*q0*vphase(x-x0,y-y0))

function make_vortex!(psi::XField,vort::PointVortex)
    @unpack xv,yv,qv = vort
    @unpack psiX,X,K,K2=psi; x,y = X
    psiX .*= @. vortexgrid(x,y',xv,yv,qv)
    return XField(psiX,X,K,K2)
end

L = 500
N = 1000

X = xvecs((L,L),(N,N))
K = kvecs((L,L),(N,N))
K2 = k2(K)
psi = one.(X[1].*X[2]') |> complex
psix = XField(psi,X,K,K2)

d = 10.0
vort1 = PointVortex(0.,d/2,1)
vort2 = PointVortex(0.,-d/2,-1)
@time make_vortex!(psix,vort1)
make_vortex!(psix,vort2)

x,y = X
psi = psix.psiX
rho = abs2.(psi)
heatmap(x,y,angle.(psi))
heatmap(x,y,rho)

vx,vy = velocity(psix)
wx = sqrt.(rho).*vx
wy = sqrt.(rho).*vy
W = wx,wy
@time Wi,Wc = helmholtz(W,psix)

ei = 0.5*abs.(Wi[1].^2 .+ Wi[2].^2)
ec = 0.5*abs.(Wc[1].^2 .+ Wc[2].^2)

x,y = X;dx = diff(x)[1]; dy = diff(y)[1]
sum(ei)*dx*dy
sum(ec)*dx*dy
heatmap(x,y,ei);xlims!(-50,50);ylims!(-50,50)
heatmap(x,y,ec);xlims!(-50,50);ylims!(-50,50)

#TODO test compressible
L = 200
N = 1000

X = xvecs((L,L),(N,N))
K = kvecs((L,L),(N,N))
K2 = k2(K)
psi = one.(X[1].*X[2]') |> complex

ktest = 2*pi
psi .+= @. 0.1*cos(ktest*X[1]*one.(X[2]'))
psix = XField(psi,X,K,K2)

rho = abs2.(psi)
heatmap(angle.(psi))
heatmap(rho)

vx,vy = velocity(psix)
wx = sqrt.(rho).*vx
wy = sqrt.(rho).*vy
@time Wi,Wc = helmholtz(wx,wy,psix)

ei = 0.5*abs.(Wi[1].^2 .+ Wi[2].^2)
ec = 0.5*abs.(Wc[1].^2 .+ Wc[2].^2)

heatmap(ei)
heatmap(ec)

#3d test
X = xvecs((1,2,3),(N,N,N))
K = kvecs((1,2,3),(N,N,N))
K2 = k2(K)
psi = randn(N,N,N) + im*randn(N,N,N)

psix = XField(psi,X,K,K2)
psik = KField(fft(psi),X,K,K2)

@time v = velocity(psix)
@time Vi,Vc = helmholtz(v,psix)


text:
using CairoMakie
using Images, FileIO
using ImageTransformations
using WindMixing
using JLD2
using WindMixing: animate_profiles_fluxes_final

u_img = rotr90(load(assetpath(joinpath(pwd(), "images_plot", "u.png"))))
v_img = rotr90(load(assetpath(joinpath(pwd(), "images_plot", "v.png"))))
T_img = rotr90(load(assetpath(joinpath(pwd(), "images_plot", "T.png"))))
uw_img = rotr90(load(assetpath(joinpath(pwd(), "images_plot", "uw.png"))))
vw_img = rotr90(load(assetpath(joinpath(pwd(), "images_plot", "vw.png"))))
wT_img = rotr90(load(assetpath(joinpath(pwd(), "images_plot", "wT.png"))))
Ri_img = rotr90(load(assetpath(joinpath(pwd(), "images_plot", "Ri.png"))))
z_img = load(assetpath(joinpath(pwd(), "images_plot", "z.png")))

axis_images = (
    u = u_img,
    v = v_img,
    T = T_img,
    z = z_img,
    uw = uw_img,
    vw = vw_img,
    wT = wT_img,
    Ri = Ri_img,
)

# EXTRACTED_FILE_PATH = joinpath(pwd(),
                        # "extracted_training_output/NDE_18sim_windcooling_windheating_18simBFGST0.8nograd_divide1f5_gradient_smallNN_leakyrelu_rate_2e-4_T0.8_1e-4_extracted.jld2")

EXTRACTED_FILE_PATH = joinpath(pwd(),
                        "extracted_training_output/NDE_3sim_diurnal_18simBFGST0.8nograd_divide1f5_gradient_smallNN_leakyrelu_rate_2e-4_T0.8_extracted.jld2")

file = jldopen(EXTRACTED_FILE_PATH)
train_files = file["training_info/train_files"]
loss_scalings = file["training_info/loss_scalings"]
train_parameters = file["training_info/parameters"]
close(file)

test_files = [
    # "wind_-5e-4_cooling_3e-8_new",
    # "wind_-2e-4_cooling_3e-8_new",
    # "wind_-2e-4_heating_-3e-8_new",
    # "wind_-5e-4_diurnal_5e-8",
    # "wind_-5e-4_diurnal_1e-8"
    # "cooling_3.5e-8_new",
    # "wind_-5e-4_new"
]

ν₀ = train_parameters["ν₀"]
ν₋ = train_parameters["ν₋"]
ΔRi = train_parameters["ΔRi"]
Riᶜ = train_parameters["Riᶜ"]
Pr = train_parameters["Pr"]

FILE_DIR = joinpath("C:\\Users\\xinle\\Documents\\OceanParameterizations.jl\\test_diurnal")
solve_oceananigans_modified_pacanowski_philander_nn(test_files, EXTRACTED_FILE_PATH, FILE_DIR,
                                                        timestep=200, convective_adjustment=false)

plot_data = NDE_profile_oceananigans(joinpath(FILE_DIR, "test_$(test_files[1])"), train_files, test_files,
                                  ν₀=ν₀, ν₋=ν₋, ΔRi=ΔRi, Riᶜ=Riᶜ, Pr=Pr,
                                  loss_scalings=loss_scalings,
                                  OUTPUT_PATH="")

# animate_profiles_fluxes_final(plot_data, axis_images,
#                     joinpath(FILE_DIR, "test_color"),
#                     animation_type="Training", n_trainings=18, training_types="Wind + Cooling, Wind + Heating", fps=60)


text:
# Maniplate entries in all metadata files to match new format.
# pwd() has to be ./2011-11-11/name/... while executing

timestamp = string(Int64(Dates.value(now())))

# Loop through all files and search for metadata specifically
for (root, dirs, files) in walkdir(".")

    meta_files = files[contains.(files, "data.txt")]
    for file in meta_files

        path = joinpath(root, file)

        # Read File
        open(path, "r+") do f
            global content = readstring(f)
        end

        # Expressions to be replaced
        regex1 = r"Micos_Positions_X_Y_Z (?<x>\d*\.?\d*) (?<y>\d*\.?\d*) (?<z>\d*\.?\d*)"

        if ismatch(regex1, content)

            # Move the file to a backup destination
            backup_dest = joinpath(root, "../meta_backup", timestamp)
            if !isdir(backup_dest)
                mkpath(backup_dest)
            end
            mv(path, joinpath(backup_dest, replace(file, "data", "backup")))

            sub1 = Base.SubstitutionString("micos_position_x\t\\g<x>\r\nmicos_position_y\t\\g<y>\r\nmicos_position_z\t\\g<z>")
            repl = replace(content, regex1, sub1)

            # Write to file
            open(path, "w") do f
                print(f, repl)
            end

            println("Changed $path")
        end
    end
end


text:
using Dates
using Pkg
using Tar

if Base.VERSION < v"1.6"
    throw(ErrorException("The `rr_capture.jl` script requires Julia 1.6 or greater"))
end

if length(ARGS) < 1
    throw(ErrorException("Usage: rr_capture.jl [command...]"))
end

const TIMEOUT = 2 * 60 * 60 # timeout in seconds

# We only use `rr` on the `tester_linux64` builder
const use_rr_if_builder_is = "tester_linux64"

const run_id = get(ENV, "BUILDKITE_JOB_ID", "unknown")
const shortcommit = get(ENV, "BUILDKITE_COMMIT", "unknown")
const builder = get(ENV, "BUILDKITE_STEP_KEY", use_rr_if_builder_is)
const use_rr = builder == use_rr_if_builder_is

@info "" run_id shortcommit builder use_rr
@info "" ARGS

# if !use_rr # TODO: uncomment this line
if true # TODO: delete this line
    @info "We will not run the tests under rr"
    p = run(`$ARGS`)
    exit(p.exitcode)
end

@info "We will run the tests under rr"

const num_cores = min(Sys.CPU_THREADS, 8, parse(Int, get(ENV, "JULIA_TEST_NUM_CORES", "8")) + 1)
@info "" num_cores

proc = nothing

new_env = copy(ENV)
mktempdir() do dir
    Pkg.activate(dir)
    Pkg.add("rr_jll")
    Pkg.add("Zstd_jll")

    rr_jll = Base.require(Base.PkgId(Base.UUID((0xe86bdf43_55f7_5ea2_9fd0_e7daa2c0f2b4)), "rr_jll"))
    zstd_jll = Base.require(Base.PkgId(Base.UUID((0x3161d3a3_bdf6_5164_811a_617609db77b4)), "Zstd_jll"))
    rr(func) = Base.invokelatest(rr_jll.rr, func; adjust_LIBPATH=false)
    rr() do rr_path
        capture_script_path = joinpath(dir, "capture_output.sh")
        loader = Sys.WORD_SIZE == 64 ? "/lib64/ld-linux-x86-64.so.2" : "/lib/ld-linux.so.2"
        open(capture_script_path, "w") do io
            write(io, """
            #!/bin/bash

            $(rr_path) record --nested=detach "\$@" > >(tee -a $(dir)/stdout.log) 2> >(tee -a $(dir)/stderr.log >&2)
            """)
        end
        chmod(capture_script_path, 0o755)

        new_env = copy(ENV)
        new_env["_RR_TRACE_DIR"] = joinpath(dir, "rr_traces")
        new_env["RR_LOG"]="all:debug"
        new_env["RR_LOG_BUFFER"]="100000"
        new_env["JULIA_RR"] = capture_script_path
        t_start = time()
        global proc = run(setenv(`$(rr_path) record --num-cores=$(num_cores) $ARGS`, new_env), (stdin, stdout, stderr); wait=false)

        # Start asynchronous timer that will kill `rr`
        @async begin
            sleep(TIMEOUT)

            # If we've exceeded the timeout and `rr` is still running, kill it.
            if isopen(proc)
                println(stderr, "\n\nProcess timed out. Signalling `rr` for force-cleanup!")
                kill(proc, Base.SIGTERM)

                # Give `rr` a chance to cleanup
                sleep(60)

                if isopen(proc)
                    println(stderr, "\n\n`rr` failed to cleanup within one minute, killing and exiting immediately!")
                    kill(proc, Base.SIGKILL)
                    exit(1)
                end
            end
        end

        # Wait for `rr` to finish, either through naturally finishing its run, or `SIGTERM`.
        # If we have to `SIGKILL`
        wait(proc)

        # On a non-zero exit code, upload the `rr` trace
        if !success(proc)
            println(stderr, "`rr` returned $(proc.exitcode), packing and uploading traces...")

            if !isdir(joinpath(dir, "rr_traces"))
                println(stderr, "No `rr_traces` directory!  Did `rr` itself fail?")
                exit(1)
            end

            # Clean up non-traces
            rm(joinpath(dir, "rr_traces", "latest-trace"))
            rm(joinpath(dir, "rr_traces", "cpu_lock"))

            # Create a directory for the pack files to go
            pack_dir = joinpath(dir, "pack")
            mkdir(pack_dir)

            # Pack all traces
            trace_dirs = [joinpath(dir, "rr_traces", f) for f in readdir(joinpath(dir, "rr_traces"))]
            filter!(isdir, trace_dirs)
            run(ignorestatus(`$(rr_path) pack --pack-dir=$pack_dir $(trace_dirs)`))

            # Tar it up
            mkpath("dumps")
            datestr = Dates.format(now(), dateformat"yyyy-mm-dd_HH_MM_SS")
            dst_path = "dumps/rr-run_$(run_id)-gitsha_$(shortcommit)-$(datestr).tar.zst"
            zstd_jll.zstdmt() do zstdp
                tarproc = open(`$zstdp -o $dst_path`, "w")
                Tar.create(dir, tarproc)
                close(tarproc.in)
            end
        end
    end
end

# Pass the exit code back up to Buildkite
if proc.termsignal != 0
    ccall(:raise, Cvoid, (Cint,), proc.termsignal)
    exit(1) # Just in case the signal did not cause an exit
else
    exit(proc.exitcode)
end

max_stars_repo_name: jonas-schulze/julia

max_issues_repo_name: jonas-schulze/julia

max_forks_repo_name: jonas-schulze/julia


text:
"Calculates counterfactuals for `d` using `impute_method`."
function getcounterfactuals(d::Data{Normalized}, impute_method::Symbol)
  cf_model = train(d, Val{impute_method}())
  cf = getoutcomes(d, cf_model)

  # Add in current prescriptions where known
  for i = 1:size(d.X, 1)
    cf[i, d.T[i]] = d.y[i]
  end

  cf
end
function getcounterfactuals(d::Data{Original}, impute_method::Symbol)
  d_norm, _, _ = normalize_data(d)
  getcounterfactuals(d_norm, impute_method)
end
function getcounterfactuals(X, y, T, impute_method::Symbol)
  getcounterfactuals(Data(X, y, T), impute_method)
end

"""
Evaluates quality of omniprescient oracle on data `d` subject to
`allowed_prescriptions`.
"""
function evaluateoracle(
    cf::Matrix{Float64},
    allowed_prescriptions::Vector{Vector{Int}}=[collect(1:size(cf, 2))
                                                for _ in 1:size(cf, 1)],
  )
  # For each sample, find the best allowed prescription
  n = length(allowed_prescriptions)
  outcomes = Vector{Float64}(undef, n)
  prescriptions = Vector{Int}(undef, n)
  for i = 1:n
    outcomes[i], idx = findmin(cf[i, allowed_prescriptions[i]])
    prescriptions[i] = allowed_prescriptions[i][idx]
  end

  outcomes, prescriptions
end

"Evaluates quality of `prescriptions` on data `d`."
function evaluateprescriptions(cf::Matrix{Float64}, prescriptions::Vector{Int})
  # For each sample, find result under the current prescription
  n = length(prescriptions)
  outcomes = Vector{Float64}(undef, n)
  for i = 1:n
    outcomes[i] = cf[i, prescriptions[i]]
  end
  outcomes
end

"Evaluates quality of baseline standard-of-care model on data `d`."
function evaluatebaseline(cf::Matrix{Float64}, d::Data)
  evaluateprescriptions(cf, d.T)
end
evaluatebaseline(cf, X, y, T) = evaluatebaseline(cf, Data(X, y, T))


text:
"""
macro gmsh(ex)

Use a `try` block to initialize `gmsh`, execute `ex`, and finalize `gmsh`
regardless of how `ex` exits.
"""
macro gmsh(ex)
    return quote
        gmsh.initialize()
        try
            set_verbosity(0)
            $(esc(ex))
        finally
            # make sure we finalize gmsh if something goes wrong
            gmsh.finalize()
        end
    end
end

"""
    domain(;dim=3)

Construct a `Domain` from the current `gmsh` model with all entities of
dimension `dim`.

!!! note
    This function assumes that `gmsh` has been initialized, and
    does not handle its finalization.
"""
function domain(;dim=3)
    Ω = Domain() # Create empty domain
    domain!(Ω;dim)
    return Ω
end

"""
    domain!(Ω::Domain;[dim=3])

Like [`domain`](@ref), but appends entities to `Ω` instead of
creating a new domain.

!!! note
    This function assumes that `gmsh` has been initialized, and does not handle its
    finalization.
"""
function domain!(Ω::Domain;dim=3)
    dim_tags = gmsh.model.getEntities(dim)
    for (_, tag) in dim_tags
        # FIXME: create a temporary GMSH_TAGS to first import all tags from
        # gmsh, then check that there are no duplicates with already existing tags?
        # if haskey(ENTITIES,(dim,tag))
        #     ent = ENTITIES[(dim,tag)]
        # else
        ent = ElementaryEntity(dim, tag)
        _fill_entity_boundary!(ent)
        # end
        push!(Ω, ent)
    end
    return Ω
end

"""
    meshgen(Ω;[dim=3])

Create a `GenericMesh` for the entities in `Ω` using the curent `gmsh` model.
Passing `d=2` will create a two-dimensional mesh by projecting the original mesh
into the `x,y` plane.

!!! danger
    This function assumes that `gmsh` has been initialized, and does not handle its
    finalization.
"""
function meshgen(Ω::Domain;dim=3)
    msh = GenericMesh{3,Float64}()
    meshgen!(msh,Ω)
    if dim == 3
        return msh
    elseif dim == 2
        return convert_to_2d(msh)
        # return msh
    else
        error("`dim` value must be `2` or `3`")
    end
end

"""
    meshgen!(msh,Ω)

Similar to [`meshgen`](@ref), but append information to `msh` instead of
creating a new mesh.

!!! danger
    This function assumes that `gmsh` has been initialized, and does not handle its
    finalization.
"""
function meshgen!(msh::GenericMesh,Ω::Domain)
    tags, coord, _ = gmsh.model.mesh.getNodes()
    # NOTE: maybe use something like "unsafe_convert" to avoid having to
    # allocate new points for the nodes?
    nodes = reinterpret(SVector{3,Float64}, coord) |> collect
    append!(msh.nodes,nodes)
    els = Mesh.elements(msh)
    e2t = ent2tags(msh)
    # Recursively populate the dictionaries
    _domain_to_mesh!(els, e2t, Ω)
    return msh
end

"""
    read_geo(fname::String;dim=3)

Read a `.geo` file and generate a [`Domain`](@ref) with all entities of
dimension `dim`.
"""
function read_geo(fname;dim=3)
    Ω = Domain() # Create empty domain
    @gmsh begin
        gmsh.open(fname)
        domain!(Ω;dim)
    end
    return Ω
end

"""
    read_msh(fname::String;dim=3)

Read `fname` and create a `Domain` and a `GenericMesh` structure with all
entities in `Ω` of dimension `dim`.
"""
function read_msh(fname;dim=3)
    Ω = Domain()
    msh = @gmsh begin
        gmsh.open(fname)
        Ω = domain(;dim)
        meshgen(Ω;dim)
    end
    return Ω,msh
end

"""
    _fill_entity_boundary!

Use the `gmsh` API to add the boundary of an `ElementaryEntity`.

This is a helper function, and should not be called by itself.
"""
function _fill_entity_boundary!(ent)
    combine  = true # FIXME: what should we use here?
    oriented = false
    dim_tags = gmsh.model.getBoundary((ent.dim, ent.tag),combine,oriented)
    for (d, t) in dim_tags
        if haskey(ENTITIES,(d,t))
            bnd = ENTITIES[(d,t)]
        else
            bnd = ElementaryEntity(d, t)
            _fill_entity_boundary!(bnd)
        end
        push!(ent.boundary, bnd)
    end
    return ent
end

"""
    _domain_to_mesh!(elements, ent2tag, Ω::Domain)

Recursively populating the dictionaries `elements` and `ent2tag`.
"""
function _domain_to_mesh!(elements, ent2tag, Ω::Domain)
    isempty(Ω) && (return elements, ent2tag)
    for ω in Ω
        _ent_to_mesh!(elements, ent2tag, ω)
    end
    Γ = skeleton(Ω)
    _domain_to_mesh!(elements, ent2tag, Γ)
end

"""
    _ent_to_mesh!(elements, ent2tag, ω::ElementaryEntity)

For each element type used to mesh `ω`:
- push into `elements::Dict` the pair `etype=>ntags`;
- push into `ent2tag::Dict` the pair `etype=>etags`;

where:
- `etype::DataType` determines the type of the element (see
    [`_type_tag_to_etype`](@ref));
- `ntags::Matrix{Int}` gives the indices of the nodes defining those
    elements;
- `etags::Vector{Int}` gives the indices of those elements in `elements`.
"""
function _ent_to_mesh!(elements, ent2tag, ω::ElementaryEntity)
    ω in keys(ent2tag) && (return elements, ent2tag)
    etypes_to_etags = OrderedDict{DataType,Vector{Int}}()
    # Loop on GMSH element types (integer)
    type_tags, _, ntagss = gmsh.model.mesh.getElements(geometric_dimension(ω),tag(ω))
    for (type_tag, ntags) in zip(type_tags, ntagss)
        _, _, _, Np, _ = gmsh.model.mesh.getElementProperties(type_tag)
        ntags = reshape(ntags, Int(Np), :)
        etype = _type_tag_to_etype(type_tag)
        if etype in keys(elements)
            etag = size(elements[etype], 2) .+ collect(1:size(ntags,2))
            ntags = hcat(elements[etype], ntags)
        else
            etag = collect(1:size(ntags, 2))
        end
        push!(elements, etype => ntags)
        push!(etypes_to_etags, etype => etag)
    end
    push!(ent2tag, ω => etypes_to_etags)
    return elements, ent2tag
end

function set_meshsize(hmax, hmin=hmax)
    gmsh.option.setNumber("Mesh.CharacteristicLengthMin", hmin)
    gmsh.option.setNumber("Mesh.CharacteristicLengthMax", hmax)
end

function set_meshorder(order)
    gmsh.option.setNumber("Mesh.ElementOrder", order)
end

function set_verbosity(i)
    gmsh.option.setNumber("General.Verbosity",i)
end

function summary(model)
    gmsh.model.setCurrent(model)
    @printf("List of entities in model `%s`: \n", model)
    @printf("|%10s|%10s|%10s|\n","name","dimension","tag")
    ents = gmsh.model.getEntities()
    # pgroups = gmsh.model.getPhysicalGroups()
    for ent in ents
        name = gmsh.model.getEntityName(ent...)
        dim, tag = ent
        @printf("|%10s|%10d|%10d|\n", name, dim, tag)
    end
    println()
end

function summary()
    models = gmsh.model.list()
    for model in models
        summary(model)
    end
end

"""
    _type_tag_to_etype(tag)

Mapping of `gmsh` element types, encoded as an integer, to the internal
equivalent of those. This function assumes `gmsh` has been initilized.
"""
function _type_tag_to_etype(tag)
    T = SVector{3,Float64} # point type
    name,dim,order,num_nodes,ref_nodes,num_primary_nodes  = gmsh.model.mesh.getElementProperties(tag)
    num_nodes = Int(num_nodes) #convert to Int64
    if occursin("Point",name)
        etype = Point3D
    elseif occursin("Line",name)
    etype = LagrangeLine{num_nodes,T}
    elseif occursin("Triangle",name)
        etype = LagrangeTriangle{num_nodes,T}
    elseif occursin("Quadrilateral",name)
        etype = LagrangeSquare{num_nodes,T}
    elseif occursin("Tetrahedron",name)
        etype = LagrangeTetrahedron{num_nodes,T}
    else
        error("gmsh element of family $name does not an internal equivalent")
    end
    return etype
end

"""
    _etype_to_type_tag(etype)

Mapping of internal element types, to the integer tag of `gmsh` elements. This
function assumes `gmsh` has been initialized.
"""
function _etype_to_type_tag(el::LagrangeElement)
    etype = typeof(el)
    tag = 1
    while true
        E   = _type_tag_to_etype(tag)
        E === etype && (return tag)
        tag = tag + 1
    end
end

"""
    struct GmshParametricEntity{M} <: AbstractEntity

An attempt to wrap a `gmsh` entity into an internal object which can access the
underlying entity parametrization.

The interface implements a somewhat *hacky* way to use `gmsh` as a `CAD` reader.
Some important points to keep in mind (which make this *hacky*):
- Each call `el(x)` entails calling a function on the `gmsh` library (linked
  dynamically). In particular the call signature for the `gmsh`'s `getValue`
  method requires passing a `Vector` for the input point, and therefore some
  unnecessary allocations are incurred. Overall, this means **this interface is inneficient**.
- Only untrimmed surfaces should be considered; thus **this interface is
  limited** to somewhat simple surfaces.
- **No guarantee is given on the quality of the parametrization**. For instance, a
  sphere may be parametrized using spherical coordinates, for which a geometric
  singularity exists at the poles. For complex surfaces, it is typically much
  better to simply mesh the surface using `gmsh` and then import the file.

!!! warning
  This is an experimental feature and should be used with caution.
"""
struct GmshParametricEntity{M} <: AbstractEntity
    tag::Int
    domain::HyperRectangle{M,Float64}
    function GmshParametricEntity{M}(tag,domain) where {M}
        dim = M
        ent = new{M}(dim,domain)
        global_add_entity!(ent)
        return ent
    # TODO: throw a warning if the surface is trimmed
    end
end

function GmshParametricEntity(dim::Int,tag::Int,model=gmsh.model.getCurrent())
    low_corner,high_corner = gmsh.model.getParametrizationBounds(dim,tag)
    rec = HyperRectangle(low_corner,high_corner)
    return GmshParametricEntity{dim}(tag,rec)
end
GmshParametricEntity(dim::Integer,tag::Integer,args...;kwargs...) = GmshParametricEntity(Int(dim),Int(tag),args...;kwargs...)

function (par::GmshParametricEntity{N})(x) where {N}
    if N === 1
        return gmsh.model.getValue(N,par.tag,x)
    elseif N===2
        return gmsh.model.getValue(N,par.tag,[x[1],x[2]])
    else
        error("got N=$N, values must be 1 or 2")
    end
end

function jacobian(psurf::GmshParametricEntity{N},s::SVector) where {N}
    if N==1
        jac = gmsh.model.getDerivative(N,psurf.tag,s)
        return SMatrix{3,1}(jac)
    elseif N==2
        jac = gmsh.model.getDerivative(N,psurf.tag,[s[1],s[2]])
        return SMatrix{3,2}(jac)
    else
        error("got N=$N, values must be 1 or 2")
    end
end

max_stars_repo_name: WaveProp/GmshSDK

max_issues_repo_name: WaveProp/GmshSDK

max_forks_repo_name: WaveProp/GmshSDK


text:
"""
Creates lossy converter model between AC and DC grid, assuming U_i is approximatley 1 numerically

```
pconv_ac[i] + pconv_dc[i] == a + b*pconv_ac
```
"""
function constraint_converter_losses(pm::AbstractDCPModel, n::Int, cnd::Int, i::Int, a, b, c, plmax)
    pconv_ac = PowerModels.var(pm, n, cnd, :pconv_ac, i)
    pconv_dc = PowerModels.var(pm, n, cnd, :pconv_dc, i)
    v = 1 #pu, assumption to approximate current
    cm_conv_ac = pconv_ac/v # can actually be negative, not a very nice model...
    if pm.setting["conv_losses_mp"] == true
        PowerModels.con(pm, n, cnd, :conv_loss)[i] = @constraint(pm.model, pconv_ac + pconv_dc == a + b*cm_conv_ac )
    else
        PowerModels.con(pm, n, cnd, :conv_loss)[i] = @constraint(pm.model, pconv_ac + pconv_dc >= a + b*cm_conv_ac )
        PowerModels.con(pm, n, cnd, :conv_loss_aux)[i] = @constraint(pm.model, pconv_ac + pconv_dc >= a - b*cm_conv_ac )
        PowerModels.con(pm, n, cnd, :conv_loss_plmax)[i] = @constraint(pm.model, pconv_ac + pconv_dc <= plmax)
    end
end
"""
Converter transformer constraints

```
p_tf_fr == -btf*(v^2)/tm*(va-vaf)
p_tf_to == -btf*(v^2)/tm*(vaf-va)
```
"""
function constraint_conv_transformer(pm::AbstractDCPModel, n::Int, cnd::Int, i::Int, rtf, xtf, acbus, tm, transformer)
    ptf_fr = PowerModels.var(pm, n, cnd, :pconv_tf_fr, i)
    ptf_to = PowerModels.var(pm, n, cnd, :pconv_tf_to, i)

    vaf = PowerModels.var(pm, n, cnd, :vaf, i)
    va = PowerModels.var(pm, n, cnd, :va, acbus)

    if transformer
        btf = imag(1/(im*xtf)) # classic DC approach to obtain susceptance form
        v = 1 # pu, assumption DC approximation
        PowerModels.con(pm, n, cnd, :conv_tf_p_fr)[i] = @constraint(pm.model, ptf_fr == -btf*(v^2)/tm*(va-vaf))
        PowerModels.con(pm, n, cnd, :conv_tf_p_to)[i] = @constraint(pm.model, ptf_to == -btf*(v^2)/tm*(vaf-va))
    else
        PowerModels.con(pm, n, cnd, :conv_tf_p_fr)[i] = @constraint(pm.model, va == vaf)
        PowerModels.con(pm, n, cnd, :conv_tf_p_to)[i] = @constraint(pm.model, ptf_fr + ptf_to  == 0)
    end
end
"""
Converter reactor constraints

```
p_pr_fr == -bc*(v^2)*(vaf-vac)
pconv_ac == -bc*(v^2)*(vac-vaf)
```
"""
function constraint_conv_reactor(pm::AbstractDCPModel, n::Int, cnd::Int, i::Int, rc, xc, reactor)
    ppr_fr = PowerModels.var(pm, n, cnd, :pconv_pr_fr, i)
    pconv_ac = PowerModels.var(pm, n, cnd, :pconv_ac, i)
    vaf = PowerModels.var(pm, n, cnd, :vaf, i)
    vac = PowerModels.var(pm, n, cnd, :vac, i)

    if reactor
        bc = imag(1/(im*xc))
        v = 1 # pu, assumption DC approximation
        PowerModels.con(pm, n, cnd, :conv_pr_p)[i] = @constraint(pm.model, ppr_fr == -bc*(v^2)*(vaf-vac))
        PowerModels.con(pm, n, cnd, :conv_pr_p_to)[i] = @constraint(pm.model, pconv_ac == -bc*(v^2)*(vac-vaf))
    else
        PowerModels.con(pm, n, cnd, :conv_pr_p)[i] =  @constraint(pm.model, vac == vaf)
        PowerModels.con(pm, n, cnd, :conv_pr_p_to)[i] = @constraint(pm.model, ppr_fr + pconv_ac  == 0)
    end
end
"""
Converter filter constraints (no active power losses)
```
p_pr_fr + p_tf_to == 0
```
"""
function constraint_conv_filter(pm::AbstractDCPModel, n::Int, cnd::Int, i::Int, bv, filter)
    ppr_fr = PowerModels.var(pm, n, cnd, :pconv_pr_fr, i)
    ptf_to = PowerModels.var(pm, n, cnd, :pconv_tf_to, i)

    PowerModels.con(pm, n, cnd, :conv_kcl_p)[i] = @constraint(pm.model,   ppr_fr + ptf_to == 0 )
end
"""
Converter current constraint (not applicable)
```
```
"""
function constraint_converter_current(pm::AbstractDCPModel, n::Int, cnd::Int, i::Int, Umax, Imax)
    # not used
end
function variable_dc_converter(pm::AbstractDCPModel; kwargs...)
    variable_converter_active_power(pm; kwargs...)
    variable_dcside_power(pm; kwargs...)
    variable_converter_filter_voltage(pm; kwargs...)
    variable_converter_internal_voltage(pm; kwargs...)
    variable_converter_to_grid_active_power(pm; kwargs...)

    variable_conv_transformer_active_power_to(pm; kwargs...)
    variable_conv_reactor_active_power_from(pm; kwargs...)
end

function variable_converter_filter_voltage(pm::AbstractDCPModel; kwargs...)
    variable_converter_filter_voltage_angle(pm; kwargs...)
end

function variable_converter_internal_voltage(pm::AbstractDCPModel; kwargs...)
    variable_converter_internal_voltage_angle(pm; kwargs...)
end
"""
Converter reactive power setpoint constraint (PF only, not applicable)
```
```
"""
function constraint_reactive_conv_setpoint(pm::AbstractDCPModel, n::Int, cnd::Int, i, qconv)
end
"""
Converter firing angle constraint (not applicable)
```
```
"""
function constraint_conv_firing_angle(pm::AbstractDCPModel, n::Int, cnd::Int, i::Int, S, P1, Q1, P2, Q2)
end
"""
Converter droop constraint (not applicable)
```
```
"""
function constraint_dc_droop_control(pm::AbstractDCPModel, n::Int, cnd::Int, i::Int, busdc_i, vref_dc, pref_dc, k_droop)
end