text:
#!/usr/bin/julia
# Renames the FASTQ files in a directory by prepending the well number based on
# a MiSeq run workbook

# Activate the proper packages
using Tk
import XLSX
using DataFrames
using Missings

# Prompt for the folder containing the FASTQs
fastq_folder = ChooseDirectory()

# Prompt for the Illumina workbook containing the well info
illumina_workbook = GetOpenFile()

# Read in the excel file
xf = XLSX.readdata(illumina_workbook, "IEM_SampleSheet!B19:D300")
fastq_ids = DataFrame(xf, :auto)
rename!(fastq_ids, Symbol.(["SampleId", "RunName", "WellNum"]))
dropmissing!(fastq_ids)
fastq_ids = string.(fastq_ids)

# Get all of the files
fastqs = readdir(fastq_folder)

# Iterate through each file
for fastq in fastqs
    # Get the full path
    fastq_path = joinpath(fastq_folder, fastq)

    # Delete if this is a JSON file
    if last(fastq,5) == ".json"
        rm(fastq_path)
        println(string("Deleting ", fastq_path))
        continue
    end

    # Pull the id based on the file name
    id = first(fastq,21)

    # Find this id in the workbook
    id_row = fastq_ids[fastq_ids.SampleId .== id,:]

    # If there are no matches, keep going
    if size(id_row)[1] < 1
        continue
    end

    # Swap the column and row
    well_num = string(last(id_row.WellNum[1],2), first(id_row.WellNum[1],1))

    # Construct the new filename
    fastq_newname = string(well_num, "_", fastq)

    # Rename the file
    mv(fastq_path, joinpath(fastq_folder, fastq_newname))
    println(string("Renaming ", fastq_path, " to ", joinpath(fastq_folder, fastq_newname)))

end

max_forks_repo_name: MillironX/wphl-bioinformatics


text:
# BinDeps integration

using BinDeps
import BinDeps: PackageManager, can_use, package_available, available_version,
    libdir, generate_steps, LibraryDependency, provider, provides, pkg_name

update_once = true

struct RPM <: PackageManager
    package
end

can_use(::Type{RPM}) = iswindows()
function package_available(p::RPM)
    global update_once
    !can_use(RPM) && return false
    pkgs = p.package
    if isa(pkgs, AbstractString)
        pkgs = [pkgs]
    end
    if (update_once::Bool)
        info("Updating WinRPM package list")
        update()
        update_once = false
    end
    return all(pkg -> (length(lookup(pkg).p) > 0), pkgs)
end

available_version(p::RPM) = lookup(p.package).p[1][xpath"version/@ver"][1]
libdir(p::RPM,dep) = joinpath(dirname(dirname(@__FILE__)), "deps", "usr", "$(Sys.ARCH)-w64-mingw32", "sys-root", "mingw", "bin")
pkg_name(p::RPM) = p.package

provider(::Type{RPM}, packages::Vector{T}; opts...) where {T<:AbstractString} = RPM(packages)
provides(::Type{RPM}, packages::AbstractArray, dep::LibraryDependency; opts...) =
    provides(provider(RPM, packages; opts...), dep; opts...)

function generate_steps(dep::LibraryDependency, h::RPM, opts)
    if get(opts, :force_rebuild, false)
        error("Will not force WinRPM to rebuild dependency \"$(dep.name)\".\n"*
              "Please make any necessary adjustments manually (This might just be a version upgrade)")
    end
    () -> install(h.package; yes=true)
end

ext: jl

max_stars_repo_name: wookay/TestWinRPM

max_issues_repo_name: wookay/TestWinRPM

max_forks_repo_name: wookay/TestWinRPM


text:
module When

macro when(cond, expr)
  if bool(eval(cond))
    eval(expr)
  end
end

end


text:
function devcount()
    ret = Ref{Cint}()
    rt.cudaGetDeviceCount(ret)
    Int(ret[])
end

function device()
    ret = Ref{Cint}()
    rt.cudaGetDevice(ret)
    Int(ret[])
end

device(dev::Integer) = (rt.cudaSetDevice(dev); dev)

# instantiate the primary contexts
const contexts = Dict{Int,CUDAdrv.CuContext}()
for dev in 0:devcount()-1
    pctx = CUDAdrv.CuPrimaryContext(CUDAdrv.CuDevice(dev))
    contexts[dev] = CUDAdrv.CuContext(pctx)
end

device_reset() = device_reset(device())

function device_reset(dev::Integer)
    # clear all items on this device from `cuda_ptrs`, so they don't get freed later
    todelete = Any[]
    for (p, pdev) in cuda_ptrs
        if pdev == dev
            finalize(p)
            push!(todelete, p)
        end
    end
    for p in todelete
        delete!(cuda_ptrs, p)
    end

    # reset the context, invalidating all derived contexts
    pctx = CUDAdrv.CuPrimaryContext(CUDAdrv.CuDevice(dev))
    CUDAdrv.unsafe_reset!(pctx, false)

    # reset the device, and instantiate a new primary context
    rt.cudaDeviceReset()
    contexts[dev] = CUDAdrv.CuContext(pctx)
end

device_synchronize() = rt.cudaDeviceSynchronize()

device_properties(dev::Integer) = (aprop = Array{rt.cudaDeviceProp}(1); rt.cudaGetDeviceProperties(aprop, dev); aprop[1])

attribute(dev::Integer, code::Integer) = (ret = Cint[0]; rt.cudaDeviceGetAttribute(ret, code, dev); Int(ret[1]))

function capability(dev::Integer)
    (attribute(dev,rt.cudaDevAttrComputeCapabilityMajor), attribute(dev,rt.cudaDevAttrComputeCapabilityMinor))
end

name(p::rt.cudaDeviceProp) = unsafe_string(convert(Ptr{UInt8}, pointer([p.name])))

# criteria = dev -> Bool
function devices(criteria::Function; nmax::Integer = typemax(Int), status=:any)
    devlist = find(map(criteria, 0:devcount().-1)).-1
    if status == :any
    elseif status == :free
        devlist = filter_free(devlist)
    end
    isempty(devlist) && error("No suitable devices found")
    devlist[1:min(nmax, length(devlist))]
end

# do syntax, f(devlist)
function devices(f::Function, criteria::Function;
                 nmax::Integer = typemax(Int), status = :any)
    devlist = devices(criteria, nmax=nmax, status=status)
    devices(f, devlist)
end

function devices(f::Function, devlist::Union{Integer,AbstractVector})
    local ret
    try
        init(devlist)
        ret = f(devlist)
    finally
        close(devlist)
    end
    ret
end

function filter_free(devlist)
    @static if libnvml != nothing
        ccall(("nvmlInit", libnvml), UInt32, ())
        freelist = Int[]
        for i in devlist
            devref = Ref{Ptr{Void}}(0)
            ccall(("nvmlDeviceGetHandleByIndex", libnvml),
                  UInt32, (UInt32, Ptr{Void}), i, devref)
            status = ccall(("nvmlDeviceGetComputeRunningProcesses", libnvml),
                           UInt32, (Ptr{Void}, UInt32, Ptr{Void}),
                           devref, Ref{UInt32}(0), C_NULL)
            if status == 0
                push!(freelist, i)
            end
        end
        ccall(("nvmlShutdown", libnvml), UInt32, ())
        return freelist
    elseif !isempty(nvidia_smi)
        smi = readstring(`$nvidia_smi`)
        if contains(smi, "No running")
            return devlist
        end
        idx = search(smi, r"[Pp]rocesses:")
        isempty(idx) && error("Neither search string found")
        lines = split(smi[last(idx):end], '\n')
        freeset = Set(devlist)
        for line in lines
            m = match(r"^\| +([0-9]+)", line)
            if m !== nothing
                length(m.captures) == 1 || error("Expected only a single capture")
                dev = parse(Int, m.captures[1])
                delete!(freeset, dev)
            end
        end
        return collect(freeset)
    else
        Base.warn_once(
        """
        Neither NVML nor nvidia-smi are available.
        Please check your setup and install either one.
        `CUDArt.filter_free` will not work.
        """)
        return devlist
    end
end

function wait_free(devlist; check_interval=10)
    f = filter_free(devlist)
    length(f) < length(devlist) && warn("Waiting for GPU devices ", devlist)
    while length(f) < length(devlist)
        sleep(check_interval)
        f = filter_free(devlist)
    end
end

# A cache of useful CUDA kernels that gets loaded and closed
# by devices(f, devlist)
immutable PtxUtils
    mod::CUDAdrv.CuModule
    fns::Dict{Any,CUDAdrv.CuFunction}
end
const global ptxdict = Dict{Integer,PtxUtils}()

function init(devlist::Union{Integer,AbstractVector})
    funcnames = ["fill_contiguous", "fill_pitched"]
    funcexts  = ["double","float","int64","uint64","int32","uint32","int16","uint16","int8","uint8"]
    datatypes = [Float64,Float32,Int64,UInt64,Int32,UInt32,Int16,UInt16,Int8,UInt8]
    utilfile  = joinpath(dirname(@__DIR__), "deps", "utils.ptx")

    # initialize all devices
    for dev in devlist
        # It has already been initialized.
        if haskey(ptxdict, dev)
            continue
        end

        device(dev)
        # allocate and destroy memory to force initialization
        free(malloc(UInt8, 1))
        # Load the utility functions.
        ptx = PtxUtils(CUDAdrv.CuModuleFile(utilfile), Dict{Any,CUDAdrv.CuFunction}())
        for func in funcnames
            for (dtype, ext) in zip(datatypes, funcexts)
                ptx.fns[(func, dtype)] = CUDAdrv.CuFunction(ptx.mod, func * "_" * ext)
            end
        end
        ptx.fns["clock_block"] = CUDAdrv.CuFunction(ptx.mod, "clock_block")
        ptxdict[dev] = ptx
    end
end

function close(devlist::Union{Integer,AbstractVector})
    for dev in devlist
        if haskey(ptxdict, dev)
            # CUDAdrv unloads the module using a finalizer
            delete!(ptxdict, dev)
        end
        device_reset(dev)
    end
end


text:
@testset "compute_nu.jl" begin

end


text:
"""
Code declaring ChunkOfSpectrum and abstract version.

Author: Eric Ford
Created: August 2020
"""

""" Abstract type for any ChunkOfSpectrum """
abstract type AbstractChunkOfSpectrum end

""" ChunkOfSpectrum for views into Spectra1DBasic or Spectra2DBasic """
struct ChunkOfSpectrum{T1<:Real,T2<:Real,T3<:Real,AA1<:AbstractArray{T1,1},AA2<:AbstractArray{T2,1},AA3<:AbstractArray{T3,1}} <: AbstractChunkOfSpectrum
    λ::AA1
    flux::AA2
    var::AA3
end

function ChunkOfSpectrum{T1,T2,T3}(λ::A1, flux::A2, var::A3) where {  T1<:Real, T2<:Real, T3<:Real, A1<:AbstractArray{T1,1}, A2<:AbstractArray{T2,1}, A3<:AbstractArray{T3,1} }
    @assert size(λ) == size(flux)
    @assert size(λ) == size(var)
    min_pixels_in_chunk = 3
    max_pixels_in_chunk = 10000
    #@assert min_pixels_in_chunk <= length(λ) <= max_pixels_in_chunk
    ChunkOfSpectrum{eltype(λ),eltype(flux),eltype(var),typeof(λ),typeof(flux),typeof(var)}(λ,flux,var)
end

function ChunkOfSpectrum(λ::A1, flux::A2, var::A3, order::Integer, pixels::AUR) where {  T1<:Real, T2<:Real, T3<:Real, A1<:AbstractArray{T1,2}, A2<:AbstractArray{T2,2}, A3<:AbstractArray{T3,2}, AUR<:AbstractUnitRange }
    @assert size(λ) == size(flux)
    @assert size(λ) == size(var)
    @assert 1 <= order <= size(λ,2)
    #println("order = ", order, "   pixels = ",pixels, "  size(λ,1) = ", size(λ,1))
    #flush(stdout)
    @assert 1 <= first(pixels) <= last(pixels) <= size(λ,1)
    ChunkOfSpectrum{T1,T2,T3}(view(λ,pixels,order),view(flux,pixels,order),view(var,pixels,order))
end

function ChunkOfSpectrum(λ::A1, flux::A2, var::A3, pixels::AUR) where {  T1<:Real, T2<:Real, T3<:Real, A1<:AbstractArray{T1,1}, A2<:AbstractArray{T2,1}, A3<:AbstractArray{T3,1}, AUR<:AbstractUnitRange }
    @assert size(λ) == size(flux)
    @assert size(λ) == size(var)
    @assert 1 <= first(pixels) <= last(pixels) <= size(λ,1)
    ChunkOfSpectrum{T1,T2,T3}(view(λ,pixels),view(flux,pixels),view(var,pixels))
end

function ChunkOfSpectrum(λ::A1, flux::A2, var::A3, loc::NamedTuple{(:pixels, :order),Tuple{AUR,I1}}) where {  T1<:Real, T2<:Real, T3<:Real, A1<:AbstractArray{T1,2}, A2<:AbstractArray{T2,2}, A3<:AbstractArray{T3,2}, AUR<:AbstractUnitRange, I1<:Integer }
    ChunkOfSpectrum(λ,flux,var,loc.order,loc.pixels)
end

function ChunkOfSpectrum(spectra::AS, order::Integer, pixels::AUR) where { AS<:AbstractSpectra2D, AUR<:AbstractUnitRange }
    ChunkOfSpectrum(spectra.λ,spectra.flux,spectra.var,order,pixels)
end

function ChunkOfSpectrum(spectra::AS, loc::NamedTuple{(:pixels, :order),Tuple{AUR,I1}}) where {  AS<:AbstractSpectra2D, AUR<:AbstractUnitRange, I1<:Integer }
    ChunkOfSpectrum(spectra.λ,spectra.flux,spectra.var,loc.order,loc.pixels)
end

function ChunkOfSpectrum(spectra::AS, pixels::AUR) where { AS<:AbstractSpectra1D, AUR<:AbstractUnitRange }
    ChunkOfSpectrum(spectra.λ,spectra.flux,spectra.var,pixels)
end

function ChunkOfSpectrum(spectra::AS, loc::NamedTuple{(:pixels, :order),Tuple{AUR,I1}}) where {  AS<:AbstractSpectra1D, AUR<:AbstractUnitRange, I1<:Integer }
    @assert loc.order == 1
    ChunkOfSpectrum(spectra.λ,spectra.flux,spectra.var,loc.pixels)
end

function empty_chunk_of_spectrum()
    ChunkOfSpectrum{Float64,Float64,Float64}(Float64[],Float64[],Float64[])
end

# TODO: Generalize in case not standard view's
get_order_index(chunk::AC) where { AC<:AbstractChunkOfSpectrum } = chunk.flux.indices[2]
get_pixels_range(chunk::AC) where { AC<:AbstractChunkOfSpectrum } = chunk.flux.indices[1]

max_stars_repo_name: RvSpectML/RvSpectMLTemplate

max_issues_repo_name: RvSpectML/RvSpectMLTemplate

max_forks_repo_name: RvSpectML/RvSpectMLTemplate


text:
vec_simple = @stochastic_model begin
    @stage 1 begin
        @parameters begin
            lb = [40.0, 20.0]
            c  = [100., 150.]
            A  = [1. 1.]
            b  = [120.]
        end
        @decision(model, x[i in 1:2] >= lb[i])
        @objective(model, Min, dot(c, x))
        @constraint(model, A * x - b in MOI.Nonpositives(1))
    end
    @stage 2 begin
        @parameters begin
            lb = [0., 0.]
            T  = [-60. 0.;
                  0.  -80.]
            W  = [6. 10.;
                  8. 5.]
        end
        @uncertain ξ[1:4]
        q  = [ξ[1], ξ[2]]
        ub = [ξ[3], ξ[4]]
        @recourse(model, lb[i] <= y[i in 1:2] <= ub[i])
        @objective(model, Max, dot(q, y))
        @constraint(model, T * x + W * y in MOI.Nonpositives(2))
    end
end

ξ₁ = @scenario ξ[1:4] = [24.0, 28.0, 500.0, 100.0] probability = 0.4
ξ₂ = @scenario ξ[1:4] = [28.0, 32.0, 300.0, 300.0] probability = 0.6

vec_simple_res = SPResult([46.67, 36.25], Dict(1 => [300., 100.], 2 => [300., 100.]), -855.83, -1518.75, 662.92, 286.92, -1445.92, -568.92)
push!(problems, (vec_simple, [ξ₁,ξ₂], vec_simple_res, "Vectorized Simple"))


text:
module AstroDates

import Dates
import Dates: year, month, day,
    hour, minute, second, millisecond,
    yearmonthday, dayofyear

const J2000 = 2.4515445e6

abstract type Calendar end

struct ProlepticJulianCalendar <: Calendar end

struct JulianCalendar <: Calendar end

struct GregorianCalendar <: Calendar end

function findyear(::ProlepticJulianCalendar, j2000day)
    -Int((-Int64(4) * j2000day - Int64(2920488)) ÷ Int64(1461))
end

function last_j2000_dayofyear(::ProlepticJulianCalendar, year)
    365 * year + (year + 1) ÷ 4 - 730123
end

isleap(::ProlepticJulianCalendar, year) = (year % 4) == 0

function findyear(::JulianCalendar, j2000day)
    -Int((-Int64(4) * j2000day - Int64(2921948)) ÷ Int64(1461))
end

function last_j2000_dayofyear(::JulianCalendar, year)
    365 * year + year ÷ 4 - 730122
end

isleap(::JulianCalendar, year) = (year % 4) == 0

function findyear(::GregorianCalendar, j2000day)
    year = Int((Int64(400) * j2000day + Int64(292194288)) ÷ Int64(146097))

    # The previous estimate is one unit too high in some rare cases
    # (240 days in the 400 years gregorian cycle, about 0.16%)
    if j2000day <= last_j2000_dayofyear(GregorianCalendar(), year - 1)
        year -= 1
    end

    year
end

function last_j2000_dayofyear(::GregorianCalendar, year)
    365 * year + year ÷ 4 - year ÷ 100 + year ÷ 400 - 730120
end

function isleap(::GregorianCalendar, year)
    year % 4 == 0 && (year % 400 == 0 || year % 100 != 0)
end

const PREVIOUS_MONTH_END_DAY_LEAP = (
    0,
    31,
    60,
    91,
    121,
    152,
    182,
    213,
    244,
    274,
    305,
    335,
)

const PREVIOUS_MONTH_END_DAY = (
    0,
    31,
    59,
    90,
    120,
    151,
    181,
    212,
    243,
    273,
    304,
    334,
)

function findmonth(dayinyear, isleap)
    offset = isleap ? 313 : 323
    dayinyear < 32 ? 1 : (10 * dayinyear + offset) ÷ 306
end

function findday(dayinyear, month, isleap)
    (!isleap && dayinyear > 365) && throw(ArgumentError("Day of year cannot be 366 for a non-leap year."))
    previous_days = isleap ? PREVIOUS_MONTH_END_DAY_LEAP : PREVIOUS_MONTH_END_DAY
    dayinyear - previous_days[month]
end

function finddayinyear(month, day, isleap)
    previous_days = isleap ? PREVIOUS_MONTH_END_DAY_LEAP : PREVIOUS_MONTH_END_DAY
    day + previous_days[month]
end

function getcalendar(year, month, day)
    calendar = GregorianCalendar()
    if year < 1583
        if year < 1
            calendar = ProlepticJulianCalendar()
        elseif year < 1582 || month < 10 || (month < 11 && day < 5)
            calendar = JulianCalendar()
        end
    end
    calendar
end

struct Date{C}
    year::Int
    month::Int
    day::Int
    Date{C}(year, month, day) where {C} = new{C::Calendar}(year, month, day)
end

@inline function Date(offset::Integer)
    calendar = GregorianCalendar()
    if offset < -152384
        if offset > -730122
            calendar = JulianCalendar()
        else
            calendar = ProlepticJulianCalendar()
        end
    end

    year = findyear(calendar, offset)
    dayinyear = offset - last_j2000_dayofyear(calendar, year - 1)

    month = findmonth(dayinyear, isleap(calendar, year))
    day = findday(dayinyear, month, isleap(calendar, year))

    Date{calendar}(year, month, day)
end

function Date(epoch::Date, offset::Int)
    Date(j2000(epoch) + offset)
end

function Date(year, month, day)
    if month < 1 || month > 12
        throw(ArgumentError("Invalid month number: $month"))
    end

    check = Date(j2000(year, month, day))
    if check.year != year || check.month != month || check.day != day
        throw(ArgumentError("Invalid date."))
    end
    Date{calendar(check)}(year, month, day)
end

function Date(year, dayinyear)
    calendar = GregorianCalendar()
    if year < 1583
        if year < 1
            calendar = ProlepticJulianCalendar()
        else
            calendar = JulianCalendar()
        end
    end
    leap = isleap(GregorianCalendar(), year)
    month = findmonth(dayinyear, leap)
    day = findday(dayinyear, month, leap)

    Date{calendar}(year, month, day)
end

year(s::Date) = s.year
month(s::Date) = s.month
day(s::Date) = s.day
calendar(s::Date{C}) where {C} = C

Base.show(io::IO, d::Date) = print(io,
                                   year(d), "-",
                                   lpad(month(d), 2, '0'), "-",
                                   lpad(day(d), 2, '0'))

Date(d::Dates.Date) = Date(Dates.year(d), Dates.month(d), Dates.day(d))
Dates.Date(d::Date) = Dates.Date(year(d), month(d), day(d))

function j2000(calendar::Calendar, year, month, day)
    last_j2000_dayofyear(calendar, year - 1) + finddayinyear(month, day, isleap(calendar, year))
end

function j2000(year, month, day)
    calendar = getcalendar(year, month, day)
    j2000(calendar, year, month, day)
end

j2000(s::Date{C}) where {C} = j2000(C, year(s), month(s), day(s))

julian(s::Date) = j2000(s) + J2000

const JULIAN_EPOCH = Date(-4712,  1,  1)
const MODIFIED_JULIAN_EPOCH = Date(1858, 11, 17)
const FIFTIES_EPOCH = Date(1950, 1, 1)
const CCSDS_EPOCH = Date(1958, 1, 1)
const GALILEO_EPOCH = Date(1999, 8, 22)
const GPS_EPOCH = Date(1980, 1, 6)
const J2000_EPOCH = Date(2000, 1, 1)
const MIN_EPOCH = Date(typemin(Int32))
const MAX_EPOCH = Date(typemax(Int32))
const UNIX_EPOCH = Date(1970, 1, 1)

struct Time{T}
    hour::Int
    minute::Int
    second::T

    function Time(hour, minute, second::T) where T
        if hour < 0 || hour > 23
            throw(ArgumentError("`hour` must be an integer between 0 and 23."))
        elseif minute < 0 || minute > 59
            throw(ArgumentError("`minute` must be an integer between 0 and 59."))
        elseif second < 0 || second >= 61.0
            throw(ArgumentError("`second` must be a float between 0 and 61."))
        end

        new{T}(hour, minute, second)
    end
end

@inline function Time(second_in_day_a, second_in_day_b)
    carry = floor(Int, second_in_day_b)
    wholeseconds = second_in_day_a + carry
    fractional = second_in_day_b - carry

    # range check
    if wholeseconds < 0 || wholeseconds > 86400
        throw(ArgumentError("Seconds are out of range"))
    end

    # extract the time components
    hour = wholeseconds ÷ 3600
    wholeseconds -= 3600 * hour
    minute = wholeseconds ÷ 60
    wholeseconds -= 60 * minute
    second = wholeseconds + fractional

    Time(hour, minute, second)
end

Time(t::Dates.Time) = Time(Dates.hour(t), Dates.minute(t),
                           Dates.second(t) + 1e-3Dates.millisecond(t))
Dates.Time(t::Time) = Dates.Time(hour(t), minute(t), second(t), millisecond(t))

const H00 = Time(0, 0, 0.0)
const H12 = Time(12, 0, 0.0)

hour(t::Time) = t.hour
minute(t::Time) = t.minute
second(::Type{Float64}, t::Time) = t.second
second(::Type{Int}, t::Time) = floor(Int, t.second)
second(t::Time) = second(Int, t)
millisecond(t::Time) = round(Int, (second(Float64, t) - second(Int, t)) * 1000)

fractionofday(t::Time) = t.second / 86400 + t.minute / 1440 + t.hour / 24

secondinday(t::Time) = t.second + 60 * t.minute + 3600 * t.hour

function Base.show(io::IO, t::Time)
    h = lpad(hour(t), 2, '0')
    m = lpad(minute(t), 2, '0')
    s = lpad(second(Int, t), 2, '0')
    f = lpad(millisecond(t), 3, '0')
    print(io, h, ":", m, ":", s, ".", f)
end

struct DateTime{C}
    date::Date{C}
    time::Time
end

Date(dt::DateTime) = dt.date
Time(dt::DateTime) = dt.time

Base.show(io::IO, dt::DateTime) = print(io, Date(dt), "T", Time(dt))

function DateTime(year, month, day, hour=0, minute=0, second=0.0)
    DateTime(Date(year, month, day), Time(hour, minute, second))
end

year(dt::DateTime) = year(Date(dt))
month(dt::DateTime) = month(Date(dt))
day(dt::DateTime) = day(Date(dt))
yearmonthday(dt::DateTime) = year(Date(dt)), month(Date(dt)), day(Date(dt))
function dayofyear(dt::DateTime{C}) where C
    leap = isleap(C, year(dt))
    finddayinyear(month(dt), day(dt), leap)
end
hour(dt::DateTime) = hour(Time(dt))
minute(dt::DateTime) = minute(Time(dt))
second(typ, dt::DateTime) = second(typ, Time(dt))
second(dt::DateTime) = second(Int, Time(dt))
millisecond(dt::DateTime) = millisecond(Time(dt))

julian(dt::DateTime) = fractionofday(Time(dt)) + julian(Date(dt))
j2000(dt::DateTime) = fractionofday(Time(dt)) + j2000(Date(dt))
julian_twopart(dt::DateTime) = julian(Date(dt)), fractionofday(Time(dt))

DateTime(dt::Dates.DateTime) = DateTime(Dates.year(dt), Dates.month(dt), Dates.day(dt),
                                        Dates.hour(dt), Dates.minute(dt),
                                        1e-3Dates.millisecond(dt) + Dates.second(dt))
function Dates.DateTime(dt::DateTime)
    y = year(dt)
    m = month(dt)
    d = day(dt)
    h = hour(dt)
    mi = minute(dt)
    s = second(Int, dt)
    ms = round((second(Float64, dt) - s) * 1000)
    Dates.DateTime(y, m, d, h, mi, s, ms)
end

end


text:
function get_initial_conditions_template(model::OperationModel)
    # This is done to avoid passing the duals but also not re-allocating the PTDF when it
    # exists
    network_model = NetworkModel(
        get_network_formulation(model.template);
        use_slacks=get_use_slacks(get_network_model(model.template)),
        PTDF=get_PTDF(get_network_model(model.template)),
    )

    ic_template = ProblemTemplate(network_model)
    for device_model in values(model.template.devices)
        base_model = get_initial_conditions_device_model(model, device_model)
        base_model.use_slacks = device_model.use_slacks
        base_model.time_series_names = device_model.time_series_names
        base_model.attributes = device_model.attributes
        set_device_model!(ic_template, base_model)
    end
    for device_model in values(model.template.branches)
        base_model = get_initial_conditions_device_model(model, device_model)
        base_model.use_slacks = device_model.use_slacks
        base_model.time_series_names = device_model.time_series_names
        base_model.attributes = device_model.attributes
        set_device_model!(ic_template, base_model)
    end
    ic_template.services = model.template.services
    return ic_template
end

function build_initial_conditions_model!(model::T) where {T <: OperationModel}
    model.internal.ic_model_container = deepcopy(get_optimization_container(model))
    ic_settings = deepcopy(model.internal.ic_model_container.settings)
    main_problem_horizon = get_horizon(ic_settings)
    # TODO: add an interface to allow user to configure initial_conditions problem
    model.internal.ic_model_container.JuMPmodel = _make_jump_model(ic_settings)
    template = get_initial_conditions_template(model)
    model.internal.ic_model_container.settings = ic_settings
    set_horizon!(ic_settings, min(INITIALIZATION_PROBLEM_HORIZON, main_problem_horizon))
    init_optimization_container!(
        model.internal.ic_model_container,
        get_network_formulation(get_template(model)),
        get_system(model),
    )
    TimerOutputs.@timeit BUILD_PROBLEMS_TIMER "Build Initialization $(get_name(model))" begin
        build_impl!(model.internal.ic_model_container, template, get_system(model))
    end
    return
end


text:
macro linklibrary_modes()
  path = normpath(Pkg.dir("FastSigmoid"),"c-src","libfastposit.so")

  esc(quote
    set_nanmode = () -> ccall( (:set_nanmode, $path), Void, (Bool,), true )
    set_infmode = () -> ccall( (:set_nanmode, $path), Void, (Bool,), false )
    set_roundstozero   = () -> ccall( (:set_underflow, $path), Void, (Bool,), true )
    set_roundsfromzero = () -> ccall( (:set_underflow, $path), Void, (Bool,), false )

    get_nanmode = () -> ccall( (:get_nanmode, $path), Bool, ())
    get_roundstozero   = () -> ccall( (:get_underflow, $path), Bool, ())
  end)
end

@linklibrary_modes

set_roundstozero()
#test that the interaction is stateful.
@test get_roundstozero();
println("testing underflow rounding mode")
#multiplication is the major operation that can cause underflow.
@linklibrary_mul(8,0)
@test pmul(0x01, 0x01) == 0x00
@test pmul(0x08, 0x08) == 0x01
@test pmul(0x08, 0x04) == 0x00
@test pmul(0x09, 0x04) == 0x01
@linklibrary_mul(8,1)
@test pmul(0x01, 0x01) == 0x00
@test pmul(0x08, 0x08) == 0x01
@test pmul(0x08, 0x04) == 0x00
@test pmul(0x09, 0x04) == 0x01
@linklibrary_mul(8,2)
@test pmul(0x01, 0x01) == 0x00
@test pmul(0x08, 0x08) == 0x01
@test pmul(0x08, 0x04) == 0x00
@test pmul(0x09, 0x04) == 0x01

@linklibrary_mul(16,0)
@test pmul(0x0001, 0x0001) == 0x0000
@test pmul(0x0080, 0x0080) == 0x0001
@test pmul(0x0080, 0x0040) == 0x0000
@test pmul(0x0081, 0x0040) == 0x0001
@linklibrary_mul(16,1)
@test pmul(0x0001, 0x0001) == 0x0000
@test pmul(0x0080, 0x0080) == 0x0001
@test pmul(0x0080, 0x0040) == 0x0000
@test pmul(0x0081, 0x0040) == 0x0001
@linklibrary_mul(16,2)
@test pmul(0x0001, 0x0001) == 0x0000
@test pmul(0x0080, 0x0080) == 0x0001
@test pmul(0x0080, 0x0040) == 0x0000
@test pmul(0x0081, 0x0040) == 0x0001

@linklibrary_mul(32,0)
@test pmul(0x0000_0001, 0x0000_0001) == 0x0000_0000
@test pmul(0x0000_8000, 0x0000_8000) == 0x0000_0001
@test pmul(0x0000_8000, 0x0000_4000) == 0x0000_0000
@test pmul(0x0000_8001, 0x0000_4000) == 0x0000_0001
@linklibrary_mul(32,1)
@test pmul(0x0000_0001, 0x0000_0001) == 0x0000_0000
@test pmul(0x0000_8000, 0x0000_8000) == 0x0000_0001
@test pmul(0x0000_8000, 0x0000_4000) == 0x0000_0000
@test pmul(0x0000_8001, 0x0000_4000) == 0x0000_0001
@linklibrary_mul(32,2)
@test pmul(0x0000_0001, 0x0000_0001) == 0x0000_0000
@test pmul(0x0000_8000, 0x0000_8000) == 0x0000_0001
@test pmul(0x0000_8000, 0x0000_4000) == 0x0000_0000
@test pmul(0x0000_8001, 0x0000_4000) == 0x0000_0001
@linklibrary_mul(32,3)
@test pmul(0x0000_0001, 0x0000_0001) == 0x0000_0000
@test pmul(0x0000_8000, 0x0000_8000) == 0x0000_0001
@test pmul(0x0000_8000, 0x0000_4000) == 0x0000_0000
@test pmul(0x0000_8001, 0x0000_4000) == 0x0000_0001

set_roundsfromzero()
#test that the interaction is stateful.
@test !get_roundstozero();
println("testing non-underflowing rounding mode")
@linklibrary_mul(8,0)
@test pmul(0x01, 0x01) == 0x01
@test pmul(0x08, 0x08) == 0x01
@test pmul(0x08, 0x04) == 0x01
@test pmul(0x09, 0x04) == 0x01
@linklibrary_mul(8,1)
@test pmul(0x01, 0x01) == 0x01
@test pmul(0x08, 0x08) == 0x01
@test pmul(0x08, 0x04) == 0x01
@test pmul(0x09, 0x04) == 0x01
@linklibrary_mul(8,2)
@test pmul(0x01, 0x01) == 0x01
@test pmul(0x08, 0x08) == 0x01
@test pmul(0x08, 0x04) == 0x01
@test pmul(0x09, 0x04) == 0x01

@linklibrary_mul(16,0)
@test pmul(0x0001, 0x0001) == 0x0001
@test pmul(0x0080, 0x0080) == 0x0001
@test pmul(0x0080, 0x0040) == 0x0001
@test pmul(0x0081, 0x0040) == 0x0001
@linklibrary_mul(16,1)
@test pmul(0x0001, 0x0001) == 0x0001
@test pmul(0x0080, 0x0080) == 0x0001
@test pmul(0x0080, 0x0040) == 0x0001
@test pmul(0x0081, 0x0040) == 0x0001
@linklibrary_mul(16,2)
@test pmul(0x0001, 0x0001) == 0x0001
@test pmul(0x0080, 0x0080) == 0x0001
@test pmul(0x0080, 0x0040) == 0x0001
@test pmul(0x0081, 0x0040) == 0x0001

@linklibrary_mul(32,0)
@test pmul(0x0000_0001, 0x0000_0001) == 0x0000_0001
@test pmul(0x0000_8000, 0x0000_8000) == 0x0000_0001
@test pmul(0x0000_8000, 0x0000_4000) == 0x0000_0001
@test pmul(0x0000_8001, 0x0000_4000) == 0x0000_0001
@linklibrary_mul(32,1)
@test pmul(0x0000_0001, 0x0000_0001) == 0x0000_0001
@test pmul(0x0000_8000, 0x0000_8000) == 0x0000_0001
@test pmul(0x0000_8000, 0x0000_4000) == 0x0000_0001
@test pmul(0x0000_8001, 0x0000_4000) == 0x0000_0001
@linklibrary_mul(32,2)
@test pmul(0x0000_0001, 0x0000_0001) == 0x0000_0001
@test pmul(0x0000_8000, 0x0000_8000) == 0x0000_0001
@test pmul(0x0000_8000, 0x0000_4000) == 0x0000_0001
@test pmul(0x0000_8001, 0x0000_4000) == 0x0000_0001
@linklibrary_mul(32,3)
@test pmul(0x0000_0001, 0x0000_0001) == 0x0000_0001
@test pmul(0x0000_8000, 0x0000_8000) == 0x0000_0001
@test pmul(0x0000_8000, 0x0000_4000) == 0x0000_0001
@test pmul(0x0000_8001, 0x0000_4000) == 0x0000_0001

################################################################################
# nan mode testing.

macro linklibrary_gt(n, es)
  path = normpath(Pkg.dir("FastSigmoid"),"c-src","libfastposit.so")
  fnname = QuoteNode(Symbol(:p, n, :e, es, :_gt))

  posittype = Symbol(:UInt, n)

  esc(:(pgt = (a, b) -> ccall( ($fnname, $path), Bool, ($posittype, $posittype), a, b )))
end
macro linklibrary_lt(n, es)
  path = normpath(Pkg.dir("FastSigmoid"),"c-src","libfastposit.so")
  fnname = QuoteNode(Symbol(:p, n, :e, es, :_gt))

  posittype = Symbol(:UInt, n)

  esc(:(plt = (a, b) -> ccall( ($fnname, $path), Bool, ($posittype, $posittype), a, b )))
end
macro linklibrary_gte(n, es)
  path = normpath(Pkg.dir("FastSigmoid"),"c-src","libfastposit.so")
  fnname = QuoteNode(Symbol(:p, n, :e, es, :_gt))

  posittype = Symbol(:UInt, n)

  esc(:(pgte = (a, b) -> ccall( ($fnname, $path), Bool, ($posittype, $posittype), a, b )))
end
macro linklibrary_lte(n, es)
  path = normpath(Pkg.dir("FastSigmoid"),"c-src","libfastposit.so")
  fnname = QuoteNode(Symbol(:p, n, :e, es, :_gt))

  posittype = Symbol(:UInt, n)

  esc(:(plte = (a, b) -> ccall( ($fnname, $path), Bool, ($posittype, $posittype), a, b )))
end
macro linklibrary_eq(n, es)
  path = normpath(Pkg.dir("FastSigmoid"),"c-src","libfastposit.so")
  fnname = QuoteNode(Symbol(:p, n, :e, es, :_eq))

  posittype = Symbol(:UInt, n)

  esc(:(peq = (a, b) -> ccall( ($fnname, $path), Bool, ($posittype, $posittype), a, b )))
end

@linklibrary_add(8,0)
@linklibrary_mul(8,0)
@linklibrary_div(8,0)
@linklibrary_gt(8,0)
@linklibrary_lt(8,0)
@linklibrary_gte(8,0)
@linklibrary_lte(8,0)
@linklibrary_eq(8,0)

set_nanmode()
@test get_nanmode()

@test padd(0x80, 0x40) == 0x80
@test padd(0x80, 0x80) == 0x80
@test pmul(0x80, 0x40) == 0x80
@test pmul(0x80, 0x00) == 0x80
@test pmul(0x80, 0x80) == 0x80
@test pdiv(0x80, 0x40) == 0x80
@test pdiv(0x80, 0x80) == 0x80
@test pdiv(0x40, 0x80) == 0x80

@test pgt(0x80, 0x40) == false
@test plt(0x80, 0x40) == false
@test pgte(0x80, 0x40) == false
@test plte(0x80, 0x40) == false
@test peq(0x80, 0x80) == false

set_infmode()
@test !get_nanmode()
@test padd(0x80, 0x40) == 0x80
#@test padd(0x80, 0x80) == 0x80  <== this should set errno.
@test pmul(0x80, 0x40) == 0x80
#@test pmul(0x80, 0x00) == 0x80 <== this should set errno.
@test pmul(0x80, 0x80) == 0x80
@test pdiv(0x80, 0x40) == 0x80
#@test pdiv(0x80, 0x80) == 0x80 <== this should set errno.
@test pdiv(0x40, 0x80) == 0x00   # note this is different.

@test pgt(0x80, 0x40) == true
@test plt(0x80, 0x40) == true
@test pgte(0x80, 0x40) == true
@test plte(0x80, 0x40) == true
@test peq(0x80, 0x80) == true


text:
a = 0; b = 1
ϵ = 0.1
nₓ = 9; nₜ = 100; t_max = 1
left_bc = PeriodicBC(); right_bc = PeriodicBC()
u₀(x) = 1.

problem = AllenCahnProblem1D(a, b, ϵ, nₓ, nₜ, t_max, left_bc, right_bc, u₀)
assemble_system!(problem)
tsm = BackwardEulerMethod()
t, u = solve(problem, tsm)

error = abs.(ones(length(u[end])) - u[end])
norm(error, Inf) < 1e-14


text:
@testset "Robots: centroidal quadruped" begin
    # TODO: add tests
    q0 = nominal_configuration(centroidal_quadruped)
    # visualizer
    vis = RoboDojo.Visualizer();
    @test visualize!(vis, box, [centroidal_quadruped], Δt=0.1);
end


text:
module Mechanisms

export getOptimizationProblem, animateFourBarMechanism,getDynamicOptimizationProblem
export PTS_VERTICAL_LINE, PTS_ELLIPTIC, PTS_PAIR_1, PTS_PAIR_2

# include("externals.jl")
# include("fourBarControl.jl")
include("testPoints.jl")
include("fourBarMechanism.jl")

end # module


text:
function decomposition_means(m_new::M, m_old::M, input_type::Symbol,
                             cond_new::Symbol, cond_old::Symbol, classes::Vector{Symbol}; forecast_string_new = "", forecast_string_old = "",
                             verbose::Symbol = :low) where M<:AbstractDSGEModel
    # Print
    println(verbose, :low, )
    info_print(verbose, :low, "Computing means of forecast decomposition...")
    println(verbose, :low, "Start time: " * string(now()))
    begin_time = time_ns()

    input_files = get_decomp_output_files(m_new, m_old, input_type, cond_new, cond_old, classes, forecast_string_new = forecast_string_new, forecast_string_old = forecast_string_old)

    for class in classes
        print(verbose, :high, "Computing " * string(class) * "...")

        # Read metadata
        variable_names = jldopen(input_files[Symbol(:decomptotal, class)]) do file
            class_long = get_class_longname(class)
            collect(keys(read(file, string(class_long) * "_indices")))
        end

        # Get to work!
        mapfcn = use_parallel_workers(m_new) ? pmap : map
        decomp_vec = mapfcn(var -> decomposition_means(m_new, m_old, input_type,
                                                       cond_new, cond_old, class, var, forecast_string_new = forecast_string_new, forecast_string_old = forecast_string_old),
                            variable_names)
        decomps = OrderedDict{Symbol, DataFrame}()
        for (var, decomp) in zip(variable_names, decomp_vec)
            decomps[var] = decomp
        end

        # Write to file
        output_file = get_decomp_mean_file(m_new, m_old, input_type, cond_new, cond_old, class, forecast_string_new = forecast_string_new, forecast_string_old = forecast_string_old)
        output_dir = dirname(output_file)
        isdir(output_dir) || mkpath(output_dir)
        jldopen(output_file, true, true, true, IOStream) do file
            write(file, "decomps", decomps)
        end
        println(verbose, :high, "wrote " * basename(output_file))
    end

    # Print
    total_mb_time     = (time_ns() - begin_time)/1e9
    total_mb_time_min = total_mb_time/60

    println(verbose, :low, "\nTotal time to compute decomposition means and bands: " * string(total_mb_time_min) * " minutes")
    println(verbose, :low, "Computation of means and bands complete: " * string(now()))
end

function decomposition_means(m_new::M, m_old::M, input_type::Symbol,
                             cond_new::Symbol, cond_old::Symbol,
                             class::Symbol, var::Symbol; forecast_string_new = "", forecast_string_old = "") where M<:AbstractDSGEModel
    # Read in dates
    input_files = get_decomp_output_files(m_new, m_old, input_type, cond_new, cond_old, [class], forecast_string_new = forecast_string_new, forecast_string_old = forecast_string_old)
    input_file = input_files[Symbol(:decomptotal, class)]
    dates = jldopen(input_file, "r") do file
        sort(collect(keys(read(file, "date_indices"))))
    end

    decomp = DataFrame(date = dates)
    for comp in [:data, :news, :shockdec, :dettrend, :para, :total]
        product = Symbol(:decomp, comp)

        input_file = input_files[Symbol(product, class)]
        #jldopen(input_file, "r") do file
        # Parse transform
        class_long = get_class_longname(class)
        transforms = load(input_file, string(class_long) * "_revtransforms")
        transform = parse_transform(transforms[var])

        # If shockdec, loop over shocks
        loopkeys = if comp == :shockdec
            shock_indices = load(input_file, "shock_indices")
            collect(keys(shock_indices))
        else
            [comp]
        end

        if comp == :shockdec
            shock_indices = load(input_file, "shock_indices")
        end

        indices = load(input_file, "$(class_long)_indices")
        var_ind = indices[var]
        for key in loopkeys
            # Read in raw output: ndraws x nperiods
            decomp_series = if comp == :shockdec
                @show key
                shock_key = shock_indices[key]
                read_forecast_series(input_file, var_ind, shock_key)
            else
                read_forecast_series(input_file, product, var_ind)
            end

            # Reverse transform
            transformed_decomp = scenario_mb_reverse_transform(decomp_series, transform, :forecast)

            # Compute mean and add to DataFrame
            decomp[key] = vec(mean(transformed_decomp, dims = 1))
        end
    end
    return decomp
end


text:
# Autogenerated wrapper script for TMatrix_jll for x86_64-apple-darwin-libgfortran3
export tmatrix_fixed_orientation, tmatrix_random_orientation

using CompilerSupportLibraries_jll
JLLWrappers.@generate_wrapper_header("TMatrix")
JLLWrappers.@declare_library_product(tmatrix_fixed_orientation, "@rpath/libtmatrixfo.dylib")
JLLWrappers.@declare_library_product(tmatrix_random_orientation, "@rpath/libtmatrixro.dylib")
function __init__()
    JLLWrappers.@generate_init_header(CompilerSupportLibraries_jll)
    JLLWrappers.@init_library_product(
        tmatrix_fixed_orientation,
        "lib/libtmatrixfo.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        tmatrix_random_orientation,
        "lib/libtmatrixro.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
function is_leap_year(year)
    # https://docs.julialang.org/en/v1/manual/control-flow/#Short-Circuit-Evaluation
    return year % 4 == 0 && (year % 100 != 0 || year % 400 == 0)
end

max_stars_repo_name: joaopalmeiro/learning-julia

max_issues_repo_name: joaopalmeiro/learning-julia

max_forks_repo_name: joaopalmeiro/learning-julia


text:
#include("../dataset_reader.jl")
const dataset2_labels = [0, 1]
const dataset2_nfeatures = 13
const dataset2_filename = "income.csv"

function get_numerical_dataset2(path::String="")
    data = get_data(joinpath(path, dataset2_filename))

    #data_percentage = 0.308 # aprox.: 10000
    #data_percentage = 0.153558 # aprox.: 5000
    data_percentage = 0.092134 # aprox.: 3000
    #data_percentage = 1.0
    data = data[1:round(Int, data_percentage*length(data))]

    workclass_dict = create_numerical_map(getindex.(data, 2))
    education_dict = create_numerical_map(getindex.(data, 3))
    marital_dict = create_numerical_map(getindex.(data, 5))
    occupation_dict = create_numerical_map(getindex.(data, 6))
    relationship_dict = create_numerical_map(getindex.(data, 7))
    race_dict = create_numerical_map(getindex.(data, 8))
    sex_dict = create_numerical_map(getindex.(data, 9))
    country_dict = create_numerical_map(getindex.(data, 13))

    conversion_pattern = [Float64, workclass_dict, education_dict, Int, marital_dict, occupation_dict, relationship_dict, race_dict, sex_dict, Int, Int, Int, country_dict, Float64]
    convert_values(data, conversion_pattern)
end

function get_dataset2(path::String=""; shuffle::Bool=false, test_perc::Float64=0.2,
                      normalise::Bool=false, separate_test_data::Bool=true, balance_classes::Bool=false)
    data = get_numerical_dataset2(path)
    separate_dataset(data, test_percentage=test_perc, shuffle=shuffle, normalise=normalise,
                     separate_test_data=separate_test_data, balance_classes=balance_classes)
end

function get_train_dataset2(path::String=""; shuffle::Bool=false, test_perc::Float64=0.2)
    data = get_numerical_dataset2(path)
    train_data, test_data = separate_data(data, test_perc; shuffle=shuffle)
    input_train_data, output_train_data = separate_input_output(train_data)
    input_train_data, output_train_data
end

function get_dataset2_setup(path::String=""; normalise::Bool=false, separate::Bool=false,
                            balance::Bool=false)
    dataset = get_dataset2(path, normalise=normalise, separate_test_data=separate, balance_classes=balance)
    if separate
        SeparatedDataset(dataset..., dataset2_labels, dataset2_nfeatures, 3)
    else
        Dataset(dataset..., dataset2_labels, dataset2_nfeatures, 3)
    end
end

#=
workclass_dict = Dict(
  "Private" => 1,
  "State-gov" => 2,
  "Without-pay" => 3,
  "Self-emp-inc" => 4,
  "Never-worked" => 5,
  "Local-gov" => 6,
  "Self-emp-not-inc" => 7,
  "Federal-gov" => 8,
  "" => 9
  ) # column 2

education_dict = Dict(
  "Preschool" => 1,
  "1st-4th" => 2,
  "5th-6th" => 3,
  "7th-8th" => 4,
  "9th" => 5,
  "10th" => 6,
  "11th" => 7,
  "12th" => 8,
  "Some-college" => 9,
  "Bachelors" => 10,
  "HS-grad" => 11,
  "Assoc-acdm" => 12,
  "Prof-school" => 13,
  "Assoc-voc" => 14,
  "Masters" => 15,
  "Doctorate" => 16
  ) # column 3

marital_dict = Dict(
  "Never-married" => 1,
  "Married-civ-spouse" => 2,
  "Married-AF-spouse" => 3,
  "Married-spouse-absent" => 4,
  "Separated" => 5,
  "Divorced" => 6,
  "Widowed" => 7
  ) # column 5

occupation_dict = Dict("" => 1,
  "Prof-specialty" => 2,
  "Farming-fishing" => 3,
  "Handlers-cleaners" => 4,
  "Transport-moving" => 5,
  "Adm-clerical" => 6,
  "Machine-op-inspct" => 7,
  "Exec-managerial" => 8,
  "Craft-repair" => 9,
  "Priv-house-serv" => 10,
  "Armed-Forces" => 11,
  "Tech-support" => 12,
  "Sales" => 13,
  "Protective-serv" => 14,
  "Other-service" => 15) # column 6

relationship_dict = Dict(
  "Wife" => 1,
  "Husband" => 2,
  "Own-child" => 3,
  "Other-relative" => 4,
  "Not-in-family" => 5,
  "Unmarried" => 6
  ) # column 7

race_dict = Dict( "Black" => 1,
 "Asian-Pac-Islander" => 2,
 "Amer-Indian-Eskimo" => 3,
 "White" => 4,
 "Other" => 5
 ) # column 8

sex_dict = Dict("Female" => 1, "Male" => 2) # column 9
=#

max_stars_repo_name: Lucas-Muniz/MNEAT

max_issues_repo_name: Lucas-Muniz/MNEAT

max_forks_repo_name: Lucas-Muniz/MNEAT


text:
using Documenter, OceanStateEstimation, Pkg
import PlutoSliderServer, CairoMakie
Pkg.precompile()

makedocs(;
    modules=[OceanStateEstimation],
    format=Documenter.HTML(),
    pages=[
        "Home" => "index.md",
        "Examples" => "examples.md",
        "Functionalities" => "API.md",
    ],
    repo="https://github.com/gaelforget/OceanStateEstimation.jl/blob/{commit}{path}#L{line}",
    sitename="OceanStateEstimation.jl",
    authors="gaelforget <gforget@mit.edu>",
    assets=String[],
)

OceanStateEstimation.CBIOMESclim_download()
OceanStateEstimation.ECCOdiags_download()
OceanStateEstimation.ECCOdiags_add("interp_coeffs")

#lst=("CBIOMES_climatogy_create.jl",)
lst=("CBIOMES_climatology_plot.jl","ECCO_standard_plots.jl")

pth_out=joinpath(@__DIR__,"build", "examples")
!isdir(pth_out) ? mkdir(pth_out) : nothing

if true #skip generating file automatically
    for i in lst
    println("PlutoSliderServer($i) started")
    fil_in=joinpath(@__DIR__,"..", "examples",i)
    fil_out=joinpath(pth_out,i[1:end-2]*"html")
    PlutoSliderServer.export_notebook(fil_in)
    mv(fil_in[1:end-2]*"html",fil_out)
    cp(fil_in,fil_out[1:end-4]*"jl")
    println("PlutoSliderServer($i) completed")
end
end

deploydocs(;
    repo="github.com/gaelforget/OceanStateEstimation.jl",
)


text:
### A Pluto.jl notebook ###
# v0.16.1

using Markdown
using InteractiveUtils

# ╔═╡ 6dd76507-be07-4c3a-9747-99eca74d3519
begin
	using Pkg
	# Pkg.develop(url = "https://github.com//likanzhan/StatisticsWithJulia.jl")
	# Pkg.activate(joinpath(Pkg.devdir(), "StatisticsWithJulia.jl"))
	Pkg.activate(Base.current_project())
	Pkg.instantiate()
	using PlutoUI, StatisticsWithJulia
	PlutoUI.TableOfContents(title = "机器学习")
end

# ╔═╡ 1ce38846-46e6-40ce-ab6d-f0040c0f9782
StatisticsWithJulia.header(title = @__FILE__, author="")

# ╔═╡ 72a095fa-c5cf-4eb1-909d-d8b49af0bdc9
StatisticsWithJulia.list_notebooks(@__FILE__)

# ╔═╡ 8a20fa4c-338c-48c9-a12a-02e6d67153df
StatisticsWithJulia.footer() #

# ╔═╡ Cell order:
# ╟─1ce38846-46e6-40ce-ab6d-f0040c0f9782
# ╟─6dd76507-be07-4c3a-9747-99eca74d3519
# ╟─72a095fa-c5cf-4eb1-909d-d8b49af0bdc9
# ╟─8a20fa4c-338c-48c9-a12a-02e6d67153df


text:
module CommunalHelperUnderwaterMusicController

using ..Ahorn, Maple

@mapdef Entity "CommunalHelper/UnderwaterMusicController" Controller(
    x::Int,
    y::Int,
    enable::Bool=false,
    dashSFX::Bool=false,
)

const placements = Ahorn.PlacementDict(
    "Underwater Music Controller (Communal Helper)" => Ahorn.EntityPlacement(
        Controller,
    ),
)

const sprite = "objects/CommunalHelper/underwaterMusicController/icon"

Ahorn.render(ctx::Ahorn.Cairo.CairoContext, entity::Controller) = Ahorn.drawSprite(ctx, sprite, 0, 0)

Ahorn.selection(entity::Controller) = Ahorn.getSpriteRectangle(sprite, Ahorn.position(entity)...)

end


text:
module _vscodeserver

include("repl.jl")
include("../debugger/debugger.jl")

function remlineinfo!(x)
    if isa(x, Expr)
        if x.head == :macrocall && x.args[2] != nothing
            id = findall(map(x -> (isa(x, Expr) && x.head == :line) || (isdefined(:LineNumberNode) && x isa LineNumberNode), x.args))
            deleteat!(x.args, id)
            for j in x.args
                remlineinfo!(j)
            end
            insert!(x.args, 2, nothing)
        else
            id = findall(map(x -> (isa(x, Expr) && x.head == :line) || (isdefined(:LineNumberNode) && x isa LineNumberNode), x.args))
            deleteat!(x.args, id)
            for j in x.args
                remlineinfo!(j)
            end
        end
    end
    x
end

using REPL, Sockets, Base64, Pkg, UUIDs
import Base: display, redisplay
import Dates
global active_module = :Main

struct InlineDisplay <: AbstractDisplay end

pid = Base.ARGS[1]

function change_module(newmodule::String, print_change = true)
    global active_module
    smods = Symbol.(split(newmodule, "."))

    val = Main
    for i = 1:length(smods)
        if isdefined(val, smods[i])
            val = getfield(val, smods[i])
        else
            println("Could not find module $newmodule")
            return
        end
    end
    expr = Meta.parse(newmodule)
    active_module = expr
    repl = Base.active_repl
    main_mode = repl.interface.modes[1]
    main_mode.prompt = string(newmodule,"> ")
    main_mode.on_done = REPL.respond(repl,main_mode; pass_empty = false) do line
        if !isempty(line)
            ex = Meta.parse(line)
            if ex isa Expr && ex.head == :module
                ret = :( Base.eval($expr, Expr(:(=), :ans, Expr(:toplevel, Meta.parse($line)))) )
            else
                ret = :( Core.eval($expr, Expr(:(=), :ans, Meta.parse($line))) )
            end
        else
            ret = :(  )
        end
        sendMsgToVscode("repl/variables", getVariables())
        return ret
    end
    print(" \r ")
    print_change && println("Changed root module to $expr")
    printstyled(string(newmodule,"> "), bold = true, color = :green)
end

function get_available_modules(m=Main, out = Module[])
    for n in names(m, all = true, imported = true)
        if isdefined(m, n) && getfield(m, n) isa Module  && !(getfield(m, n) in out)
            M = getfield(m, n)
            push!(out, M)
            get_available_modules(M, out)
        end
    end
    out
end

function getVariables()
    M = @__MODULE__
    variables = []
    msg = ""
    for n in names(M)
        !isdefined(M, n) && continue
        x = getfield(M, n)
        t = typeof(x)
        msg = string(msg, ";", n, "::", t)
    end
    return msg
end

function generate_pipe_name(name)
    if Sys.iswindows()
        "\\\\.\\pipe\\vscode-language-julia-$name-$pid"
    elseif Sys.isunix()
        joinpath(tempdir(), "vscode-language-julia-$name-$pid")
    end
end

!(Sys.isunix() || Sys.iswindows()) && error("Unknown operating system.")

pipename_fromrepl = generate_pipe_name("fromrepl")
pipename_torepl = generate_pipe_name("torepl")

if issocket(pipename_torepl)
    rm(pipename_torepl)
end

function sendMsgToVscode(cmd, payload)
    println(conn, cmd, ":", sizeof(payload))
    print(conn, payload)
end

function ends_with_semicolon(x)
    return REPL.ends_with_semicolon(split(x,'\n',keepempty = false)[end])
end

@async begin
    server = listen(pipename_torepl)
    global conn = connect(pipename_fromrepl)
    while true
        sock = accept(server)
        header = readline(sock)
        cmd, payload_size_asstring = split(header, ':')
        payload_size = parse(Int, payload_size_asstring)
        payload = read(sock, payload_size)
        if cmd == "repl/include"
            text = String(payload)
            cmod = Core.eval(active_module)
            ex = Expr(:call, :include, strip(text, '\n'))
            cmod.eval(ex)
        elseif cmd == "repl/getVariables"
            sendMsgToVscode("repl/variables", getVariables())
        elseif cmd == "debug/info"
            @info "RECEIVED A debug/info message"
            @info "With payload_size=$payload_size"
            @info String(payload)
        elseif cmd == "repl/runcode"
            payload_as_string = String(payload)
            end_first_line_pos = findfirst("\n", payload_as_string)[1]
            end_second_line_pos = findnext("\n", payload_as_string, end_first_line_pos+1)[1]

            source_filename = payload_as_string[1:end_first_line_pos-1]
            code_line, code_column = parse.(Int, split(payload_as_string[end_first_line_pos+1:end_second_line_pos-1], ':'))
            source_code = payload_as_string[end_second_line_pos+1:end]

            hideprompt() do
                if isdefined(Main, :Revise) && isdefined(Main.Revise, :revise) && Main.Revise.revise isa Function
                    let mode = get(ENV, "JULIA_REVISE", "auto")
                        mode == "auto" && Main.Revise.revise()
                    end
                end
                for (i,line) in enumerate(eachline(IOBuffer(source_code)))
                    if i==1
                        printstyled("julia> ", color=:green)
                        print(' '^code_column)
                    else
                        # Indent by 7 so that it aligns with the julia> prompt
                        print(' '^7)
                    end

                    println(line)
                end

                try
                    withpath(source_filename) do
                        res = Base.invokelatest(include_string, Main, '\n'^code_line * ' '^code_column *  source_code, source_filename)

                        if res !== nothing && !ends_with_semicolon(source_code)
                            Base.invokelatest(display, res)
                        end
                    end
                catch err
                    Base.display_error(stderr, err, catch_backtrace())
                end
            end
        elseif cmd == "repl/startdebugger"
            hideprompt() do
                payload_as_string = String(payload)
                try
                    VSCodeDebugger.startdebug(payload_as_string)
                catch err
                    Base.display_error(err, catch_backtrace())
                end
            end
        end
    end
end

function display(d::InlineDisplay, ::MIME{Symbol("image/png")}, x)
    payload = stringmime(MIME("image/png"), x)
    sendMsgToVscode("image/png", payload)
end

displayable(d::InlineDisplay, ::MIME{Symbol("image/png")}) = true

function display(d::InlineDisplay, ::MIME{Symbol("image/svg+xml")}, x)
    payload = stringmime(MIME("image/svg+xml"), x)
    sendMsgToVscode("image/svg+xml", payload)
end

displayable(d::InlineDisplay, ::MIME{Symbol("image/svg+xml")}) = true

function display(d::InlineDisplay, ::MIME{Symbol("text/html")}, x)
    payload = stringmime(MIME("text/html"), x)
    sendMsgToVscode("text/html", payload)
end

displayable(d::InlineDisplay, ::MIME{Symbol("text/html")}) = true

function display(d::InlineDisplay, ::MIME{Symbol("juliavscode/html")}, x)
    payload = stringmime(MIME("juliavscode/html"), x)
    sendMsgToVscode("juliavscode/html", payload)
end

Base.Multimedia.istextmime(::MIME{Symbol("juliavscode/html")}) = true

displayable(d::InlineDisplay, ::MIME{Symbol("juliavscode/html")}) = true

function display(d::InlineDisplay, ::MIME{Symbol("application/vnd.vegalite.v2+json")}, x)
    payload = stringmime(MIME("application/vnd.vegalite.v2+json"), x)
    sendMsgToVscode("application/vnd.vegalite.v2+json", payload)
end

displayable(d::InlineDisplay, ::MIME{Symbol("application/vnd.vegalite.v2+json")}) = true

function display(d::InlineDisplay, ::MIME{Symbol("application/vnd.vegalite.v3+json")}, x)
    payload = stringmime(MIME("application/vnd.vegalite.v3+json"), x)
    sendMsgToVscode("application/vnd.vegalite.v3+json", payload)
end

displayable(d::InlineDisplay, ::MIME{Symbol("application/vnd.vegalite.v3+json")}) = true

function display(d::InlineDisplay, ::MIME{Symbol("application/vnd.vegalite.v4+json")}, x)
    payload = stringmime(MIME("application/vnd.vegalite.v4+json"), x)
    sendMsgToVscode("application/vnd.vegalite.v4+json", payload)
end

displayable(d::InlineDisplay, ::MIME{Symbol("application/vnd.vegalite.v4+json")}) = true

function display(d::InlineDisplay, ::MIME{Symbol("application/vnd.vega.v3+json")}, x)
    payload = stringmime(MIME("application/vnd.vega.v3+json"), x)
    sendMsgToVscode("application/vnd.vega.v3+json", payload)
end

displayable(d::InlineDisplay, ::MIME{Symbol("application/vnd.vega.v3+json")}) = true

function display(d::InlineDisplay, ::MIME{Symbol("application/vnd.vega.v4+json")}, x)
    payload = stringmime(MIME("application/vnd.vega.v4+json"), x)
    sendMsgToVscode("application/vnd.vega.v4+json", payload)
end

displayable(d::InlineDisplay, ::MIME{Symbol("application/vnd.vega.v4+json")}) = true

function display(d::InlineDisplay, ::MIME{Symbol("application/vnd.vega.v5+json")}, x)
    payload = stringmime(MIME("application/vnd.vega.v5+json"), x)
    sendMsgToVscode("application/vnd.vega.v5+json", payload)
end

displayable(d::InlineDisplay, ::MIME{Symbol("application/vnd.vega.v5+json")}) = true

function display(d::InlineDisplay, ::MIME{Symbol("application/vnd.plotly.v1+json")}, x)
    payload = stringmime(MIME("application/vnd.plotly.v1+json"), x)
    sendMsgToVscode("application/vnd.plotly.v1+json", payload)
end

displayable(d::InlineDisplay, ::MIME{Symbol("application/vnd.dataresource+json")}) = true

function display(d::InlineDisplay, ::MIME{Symbol("application/vnd.dataresource+json")}, x)
    payload = stringmime(MIME("application/vnd.dataresource+json"), x)
    sendMsgToVscode("application/vnd.dataresource+json", payload)
end

Base.Multimedia.istextmime(::MIME{Symbol("application/vnd.dataresource+json")}) = true

displayable(d::InlineDisplay, ::MIME{Symbol("application/vnd.plotly.v1+json")}) = true

function display(d::InlineDisplay, x)
    if showable("application/vnd.vegalite.v4+json", x)
        display(d,"application/vnd.vegalite.v4+json", x)
    elseif showable("application/vnd.vegalite.v3+json", x)
        display(d,"application/vnd.vegalite.v3+json", x)
    elseif showable("application/vnd.vegalite.v2+json", x)
        display(d,"application/vnd.vegalite.v2+json", x)
    elseif showable("application/vnd.vega.v5+json", x)
        display(d,"application/vnd.vega.v5+json", x)
    elseif showable("application/vnd.vega.v4+json", x)
        display(d,"application/vnd.vega.v4+json", x)
    elseif showable("application/vnd.vega.v3+json", x)
        display(d,"application/vnd.vega.v3+json", x)
    elseif showable("application/vnd.plotly.v1+json", x)
        display(d,"application/vnd.plotly.v1+json", x)
    elseif showable("juliavscode/html", x)
        display(d,"juliavscode/html", x)
    # elseif showable("text/html", x)
    #     display(d,"text/html", x)
    elseif showable("image/svg+xml", x)
        display(d,"image/svg+xml", x)
    elseif showable("image/png", x)
        display(d,"image/png", x)
    else
        throw(MethodError(display,(d,x)))
    end
end

function _display(d::InlineDisplay, x)
    if showable("application/vnd.dataresource+json", x)
        display(d, "application/vnd.dataresource+json", x)
    else
        display(d, x)
    end
end

if length(Base.ARGS) >= 3 && Base.ARGS[3] == "true"
    atreplinit(i->Base.Multimedia.pushdisplay(InlineDisplay()))
end

# Load revise?
load_revise = Base.ARGS[2] == "true"

const tabletraits_uuid = UUIDs.UUID("3783bdb8-4a98-5b6b-af9a-565f29a5fe9c")
const datavalues_uuid = UUIDs.UUID("e7dc6d0d-1eca-5fa6-8ad6-5aecde8b7ea5")

global _isiterabletable = i -> false
global _getiterator = i -> i

function pkgload(pkg)
    if pkg.uuid==tabletraits_uuid
        x = Base.require(pkg)

        global _isiterabletable = x.isiterabletable
        global _getiterator = x.getiterator
    elseif pkg.uuid==datavalues_uuid
        x = Base.require(pkg)

        eval(quote
            function JSON_print_escaped(io, val::$(x.DataValue))
                $(x.isna)(val) ? print(io, "null") : JSON_print_escaped(io, val[])
            end

            julia_type_to_schema_type(::Type{T}) where {S, T<:$(x.DataValue){S}} = julia_type_to_schema_type(S)
        end)
    end
end

push!(Base.package_callbacks, pkgload)

struct CachedDataResourceString
    content::String
end
Base.show(io::IO, ::MIME"application/vnd.dataresource+json", source::CachedDataResourceString) = print(io, source.content)
Base.showable(::MIME"application/vnd.dataresource+json", dt::CachedDataResourceString) = true

function JSON_print_escaped(io, val::AbstractString)
    print(io, '"')
    for c in val
        if c=='"' || c=='\\'
            print(io, '\\')
            print(io, c)
        elseif c=='\b'
            print(io, '\\')
            print(io, 'b')
        elseif c=='\f'
            print(io, '\\')
            print(io, 'f')
        elseif c=='\n'
            print(io, '\\')
            print(io, 'n')
        elseif c=='\r'
            print(io, '\\')
            print(io, 'r')
        elseif c=='\t'
            print(io, '\\')
            print(io, 't')
        else
            print(io, c)
        end
    end
    print(io, '"')
end

function JSON_print_escaped(io, val)
    print(io, '"')
    print(io, val)
    print(io, '"')
end

function JSON_print_escaped(io, val::Missing)
    print(io, "null")
end

julia_type_to_schema_type(::Type{T}) where {T} = "string"
julia_type_to_schema_type(::Type{T}) where {T<:AbstractFloat} = "number"
julia_type_to_schema_type(::Type{T}) where {T<:Integer} = "integer"
julia_type_to_schema_type(::Type{T}) where {T<:Bool} = "boolean"
julia_type_to_schema_type(::Type{T}) where {T<:Dates.Time} = "time"
julia_type_to_schema_type(::Type{T}) where {T<:Dates.Date} = "date"
julia_type_to_schema_type(::Type{T}) where {T<:Dates.DateTime} = "datetime"
julia_type_to_schema_type(::Type{T}) where {T<:AbstractString} = "string"

function printdataresource(io::IO, source)
    if Base.IteratorEltype(source) isa Base.EltypeUnknown
        first_el = first(source)
        col_names = String.(propertynames(first_el))
        col_types = [fieldtype(typeof(first_el), i) for i=1:length(col_names)]
    else
        col_names = String.(fieldnames(eltype(source)))
        col_types = [fieldtype(eltype(source), i) for i=1:length(col_names)]
    end

    print(io, "{")

    JSON_print_escaped(io, "schema")
    print(io, ": {")
    JSON_print_escaped(io, "fields")
    print(io, ":[")
    for i=1:length(col_names)
        if i>1
            print(io, ",")
        end

        print(io, "{")
        JSON_print_escaped(io, "name")
        print(io, ":")
        JSON_print_escaped(io, col_names[i])
        print(io, ",")
        JSON_print_escaped(io, "type")
        print(io, ":")
        JSON_print_escaped(io, julia_type_to_schema_type(col_types[i]))
        print(io, "}")
    end
    print(io, "]},")

    JSON_print_escaped(io, "data")
    print(io, ":[")

    for (row_i, row) in enumerate(source)
        if row_i>1
            print(io, ",")
        end

        print(io, "{")
        for col in 1:length(col_names)
            if col>1
                print(io, ",")
            end
            JSON_print_escaped(io, col_names[col])
            print(io, ":")
            # TODO This is not type stable, should really unroll the loop in a generated function
            JSON_print_escaped(io, row[col])
        end
        print(io, "}")
    end

    print(io, "]}")
end

function remove_lln!(ex::Expr)
    for i in length(ex.args):-1:1
        if ex.args[i] isa LineNumberNode
            deleteat!(ex.args, i)
        elseif ex.args[i] isa Expr
            remove_lln!(ex.args[i])
        end
    end
end

end

function vscodedisplay(x)
    if showable("application/vnd.dataresource+json", x)
        _vscodeserver._display(_vscodeserver.InlineDisplay(), x)
    elseif _vscodeserver._isiterabletable(x)===true
        buffer = IOBuffer()
        io = IOContext(buffer, :compact=>true)
        _vscodeserver.printdataresource(io, _vscodeserver._getiterator(x))
        buffer_asstring = _vscodeserver.CachedDataResourceString(String(take!(buffer)))
        _vscodeserver._display(_vscodeserver.InlineDisplay(), buffer_asstring)
    elseif _vscodeserver._isiterabletable(x)===missing
        try
            buffer = IOBuffer()
            io = IOContext(buffer, :compact=>true)
            _vscodeserver.printdataresource(io, _vscodeserver._getiterator(x))
            buffer_asstring = _vscodeserver.CachedDataResourceString(String(take!(buffer)))
            _vscodeserver._display(_vscodeserver.InlineDisplay(), buffer_asstring)
        catch err
            _vscodeserver._display(_vscodeserver.InlineDisplay(), x)
        end
    else
        _vscodeserver._display(_vscodeserver.InlineDisplay(), x)
    end
end

vscodedisplay() = i -> vscodedisplay(i)

if _vscodeserver.load_revise
    try
        @eval using Revise
        Revise.async_steal_repl_backend()
    catch err
    end
end

macro enter(command)
    _vscodeserver.remove_lln!(command)
    :(_vscodeserver.sendMsgToVscode("debugger/enter", $(string(command))))
end

macro run(command)
    _vscodeserver.remove_lln!(command)
    :(_vscodeserver.sendMsgToVscode("debugger/run", $(string(command))))
end

max_stars_repo_name: oschulz/julia-vscode

max_issues_repo_name: oschulz/julia-vscode

max_forks_repo_name: oschulz/julia-vscode


text:
using ExecutableSpecifications
using ExecutableSpecifications:
    StepDefinitionMatcher,
    Executor, executescenario, @expect, SuccessfulStepExecution,
    issuccess, beforescenario, afterscenario, FromSourceExecutionEnvironment

mutable struct FakeExecutionEnvironment <: ExecutableSpecifications.ExecutionEnvironment
    afterscenariowasexecuted::Bool

    FakeExecutionEnvironment() = new(false)
end

function ExecutableSpecifications.beforescenario(::FakeExecutionEnvironment, context::StepDefinitionContext, scenario::Gherkin.Scenario)
    context[:beforescenariowasexecuted] = true
end

function ExecutableSpecifications.afterscenario(f::FakeExecutionEnvironment, context::StepDefinitionContext, scenario::Gherkin.Scenario)
    context[:afterscenariowasexecuted] = true
    f.afterscenariowasexecuted = true
end

struct SingleStepDefinitionMatcher <: StepDefinitionMatcher
    stepbody::Function
end

function ExecutableSpecifications.findstepdefinition(
        s::SingleStepDefinitionMatcher,
        step::ExecutableSpecifications.Gherkin.ScenarioStep)
    stepdefinition = context -> begin
        s.stepbody(context)
        SuccessfulStepExecution()
    end
    StepDefinition("some text", stepdefinition, StepDefinitionLocation("", 0))
end

@testset "Execution Environment" begin
    @testset "Scenario Execution Environment" begin
        @testset "beforescenario is defined; beforescenario is executed before the scenario" begin
            # Arrange
            env = FakeExecutionEnvironment()

            # This step definition tests that the symbol :beforescenariowasexecuted is present in
            # the context at execution time. This symbol should be set by the
            # FakeExecutionEnvironment.
            stepdefmatcher = SingleStepDefinitionMatcher(context -> @assert context[:beforescenariowasexecuted])
            executor = Executor(stepdefmatcher; executionenv=env)
            scenario = Scenario("Description", [], [Given("")])

            # Act
            result = executescenario(executor, scenario)

            # Assert
            @test issuccess(result.steps[1])
        end

        @testset "beforescenario is the default noop; beforescenario is not executed before the scenario" begin
            # Arrange
            stepdefmatcher = SingleStepDefinitionMatcher(context -> @assert !haskey(context, :beforescenariowasexecuted))
            executor = Executor(stepdefmatcher)
            scenario = Scenario("Description", [], [Given("")])

            # Act
            result = executescenario(executor, scenario)

            # Assert
            @test issuccess(result.steps[1])
        end

        @testset "afterscenario is defined; afterscenario has not been called at scenario execution" begin
            # Arrange
            env = FakeExecutionEnvironment()

            # This step definition tests that the symbol :beforescenariowasexecuted is present in
            # the context at execution time. This symbol should be set by the
            # FakeExecutionEnvironment.
            stepdefmatcher = SingleStepDefinitionMatcher(context -> @assert !haskey(context, :afterscenariowasexecuted))
            executor = Executor(stepdefmatcher; executionenv=env)
            scenario = Scenario("Description", [], [Given("")])

            # Act
            result = executescenario(executor, scenario)

            # Assert
            @test issuccess(result.steps[1])
        end

        @testset "afterscenario is defined; afterscenario is called" begin
            # Arrange
            env = FakeExecutionEnvironment()

            # This step definition tests that the symbol :beforescenariowasexecuted is present in
            # the context at execution time. This symbol should be set by the
            # FakeExecutionEnvironment.
            stepdefmatcher = SingleStepDefinitionMatcher(context -> nothing)
            executor = Executor(stepdefmatcher; executionenv=env)
            scenario = Scenario("Description", [], [Given("")])

            # Act
            result = executescenario(executor, scenario)

            # Assert
            @test env.afterscenariowasexecuted
        end
    end

    @testset "FromSourceExecutionEnvironment" begin
        @testset "beginscenario is defined in source; beginscenario is executed" begin
            # Arrange
            env = FromSourceExecutionEnvironment("""
                using ExecutableSpecifications

                @beforescenario begin
                    context[:beforescenariowasexecuted] = true
                end
            """)

            context = StepDefinitionContext()
            scenario = Gherkin.Scenario("", [], [])

            # Act
            beforescenario(env, context, scenario)

            # Assert
            @test context[:beforescenariowasexecuted]
        end

        @testset "No beginscenario is defined in source; beginscenario is a noop" begin
            # Arrange
            env = FromSourceExecutionEnvironment("")

            context = StepDefinitionContext()
            scenario = Gherkin.Scenario("", [], [])

            # Act
            beforescenario(env, context, scenario)

            # Assert
            @test !haskey(context, :beforescenariowasexecuted)
        end

        @testset "afterscenario is defined in source; afterscenario is executed" begin
            # Arrange
            env = FromSourceExecutionEnvironment("""
                using ExecutableSpecifications

                @afterscenario begin
                    context[:afterscenariowasexecuted] = true
                end
            """)

            context = StepDefinitionContext()
            scenario = Gherkin.Scenario("", [], [])

            # Act
            afterscenario(env, context, scenario)

            # Assert
            @test context[:afterscenariowasexecuted]
        end

        @testset "No afterscenario is defined in source; afterscenario is a noop" begin
            # Arrange
            env = FromSourceExecutionEnvironment("")

            context = StepDefinitionContext()
            scenario = Gherkin.Scenario("", [], [])

            # Act
            afterscenario(env, context, scenario)

            # Assert
            @test !haskey(context, :afterscenariowasexecuted)
        end
    end
end


text:
@testset "Run" begin
    @testset "constructors" begin
        tuning_run = Run()
        @test typeof(tuning_run) <: Run
    end
end


text:
# This file is part of the IntervalArithmetic.jl package; MIT licensed

"""
    widened(x)

widens the lowest and highest bounds of `x` to the previous and next representable floating-point numbers, respectively.

It implements `nextOut` from the IEEE 1788-2015 standard.
"""
widened(x::Interval{T}) where {T<:AbstractFloat} =
    unsafe_interval(prevfloat(x.lo), nextfloat(x.hi))
widened(x::EmptyInterval) = x

# unsafe method, only for use with checked input
signed_infinity(x::T) where {T<:AbstractFloat} = sign(x) * T(Inf) |>
    x -> unsafe_interval(prevfloat(x), nextfloat(x))

# common thin intervals
zero(a::I) where {T, I <: ValidInterval{T}} = unsafe_interval(zero(T))
zero(::Type{I}) where {T, I <: ValidInterval{T}} = unsafe_interval(zero(T))
one(a::I) where {T, I <: ValidInterval{T}} = unsafe_interval(one(T))
one(::Type{I}) where {T, I <: ValidInterval{T}} = unsafe_interval(one(T))
#eps(a::Interval) = unsafe_interval(max(eps(a.lo), eps(a.hi)))
eps(::Type{ValidInterval{T}}) where T = unsafe_interval(eps(T))

copysign(::ValidInterval{T}, ::ValidInterval{T}) where T = emptyinterval(T)
copysign(a::Interval{T}, b::Interval{T}) where T = abs(a) * sign(b)

"""
    valtype(::AbstractInterval{T}) where T = T

Returns the type of the bounds of the interval.

### Example

```julia
julia> valtype(Interval(1.0, 2.0))
Float64
julia> valtype(emptyinterval(Float32))
Float32
```
"""
valtype(::AbstractInterval{T}) where T = T

"""
Computes the integer hash code for an `Interval` using the method for composite types used in `AutoHashEquals.jl`
"""
hash(x::Interval, h::UInt) = hash(x.hi, hash(x.lo, hash(:Interval, h)))
hash(::EmptyInterval, h::UInt) = hash(:EmptyInterval, h)
hash(::NotAnInterval, h::UInt) = hash(:NotAnInterval, h)


text:
# Generate a smooth function whose roots are the solutions of the MCP.
# More precisely, this function is:
#
# x -> phiminus(phiplus(f(x), x-upper), x-lower)
#
# where
#  phiplus(u, v) = u + v + sqrt(u^2+v^2)
#  phiminus(u, v) = u + v - sqrt(u^2+v^2)
# are the Fischer functions.
#
# Note that Miranda and Fackler use the opposite sign convention for the MCP,
# hence the difference in the smooth function.
struct MCP
end
function mcp_smooth(df::OnceDifferentiable,
                    lower::Vector, upper::Vector)

    function f!(F, x)
        value!!(df, F, x)
        for i = 1:length(x)
            if  isfinite.(upper[i])
                F[i] += (x[i]-upper[i]) + sqrt(F[i]^2+(x[i]-upper[i])^2)
            end
            if  isfinite.(lower[i])
                F[i] += (x[i]-lower[i]) - sqrt(F[i]^2+(x[i]-lower[i])^2)
            end
        end
    end

    function j!(J, x)
        F = copy(x)
        value_jacobian!!(df, F, J, x)

        # Derivatives of phiplus
        sqplus = sqrt.(F.^2 .+ (x .- upper).^2)

        dplus_du = 1 .+ F./sqplus

        dplus_dv = copy(x)
        for i = 1:length(x)
            if isfinite.(upper[i])
                dplus_dv[i] = 1 + (x[i]-upper[i])/sqplus[i]
            else
                dplus_dv[i] = 0
            end
        end

        # Derivatives of phiminus
        phiplus = copy(F)
        for i = 1:length(x)
            if isfinite(upper[i])
                phiplus[i] += (x[i]-upper[i]) + sqplus[i]
            end
        end

        sqminus = sqrt.(phiplus.^2 .+ (x .- lower).^2)

        dminus_du = 1 .- phiplus ./ sqminus

        dminus_dv = copy(x)
        for i = 1:length(x)
            if isfinite.(lower[i])
                dminus_dv[i] = 1 - (x[i]-lower[i])/sqminus[i]
            else
                dminus_dv[i] = 0
            end
        end

        # Final computations
        for i = 1:length(x)
            for j = 1:length(x)
                J[i,j] *= dminus_du[i]*dplus_du[i]
            end
            J[i,i] += dminus_dv[i] + dminus_du[i]*dplus_dv[i]
        end
    end
    return OnceDifferentiable(f!, j!, copy(df.x_f), copy(df.F))
end

# Generate a function whose roots are the solutions of the MCP.
# More precisely, this function is:
#
# x -> min(max(f(x), x-upper), x-lower)
#
# Note that Miranda and Fackler use the opposite sign convention for the MCP,
# hence the difference in the function.
function mcp_minmax(df::OnceDifferentiable,
                    lower::Vector, upper::Vector)
    function f!(F, x)
        value!!(df, F, x)
        for i = 1:length(x)
            if F[i] < x[i]-upper[i]
                F[i] = x[i]-upper[i]
            end
            if F[i] > x[i]-lower[i]
                F[i] = x[i]-lower[i]
            end
        end
    end

    function j!(J, x)
        F = copy(x)
        value_jacobian!!(df, F, J, x)
        for i = 1:length(x)
            if F[i] < x[i]-upper[i] || F[i] > x[i]-lower[i]
                for j = 1:length(x)
                    J[i,j] = (i == j ? 1.0 : 0.0)
                end
            end
        end
    end
    return OnceDifferentiable(f!, j!, copy(df.x_f), copy(df.F))
end


text:
ENV["MATRIXDEPOT_URL_REDIRECT"] = length(VERSION.prerelease) == 0 ? "1" : "1"
const basedir = tempname() # tests start with a clean data directory
ENV["MATRIXDEPOT_DATA"] = abspath(basedir, "data")

using MatrixDepot
using Test
using LinearAlgebra
using SparseArrays

import MatrixDepot: DataError

# that will download the index files if necessary and initialize internal data
MatrixDepot.init()

@testset "MatrixDepot.jl" begin

include("generators.jl")
include("include_generator.jl")

@testset "MatrixDepot simulate remote matrix tests" begin
    tests = [
            "download",
            "common",
            "number",
            "property",
            "downloadmm",
            ]

    @testset "$t" for t in tests
        tp = joinpath(@__DIR__(), "$(t).jl")
        println("running $(tp) ...")
        include(tp)
        println("finished $(tp)")
    end
end

@testset "MatrixDepot real remote matrix tests" begin
    println("running remote.jl ...")
    include("remote.jl")
    println("finished remote.jl")
end

end # testset MatrixDepot.jl


text:
module sb

export OhmicSBParams, SBModel,
       EhrenfestOps, RunEhrenfest!,
       ECEIDOps, RunECEID!,
       CEIDOps, RunCEID!,
       H, F, K, comm, acomm, IOHBAR

include("physical_constants.jl")

include("matrices.jl")

include("ohmicsb.jl")

include("ehrenfest.jl")

include("ECEID.jl")

include("CEID.jl")

end # module


text:
# This file is part of RLinearAlgebra.jl

module ProceduralTestRLSSolver

using Test, RLinearAlgebra, Random

@testset "Randomized Linear Solver -- Procedural" begin

    # Check if solver encapsulation is defined
    @test @isdefined RLSSolver

    # Verify required fields
    @test :sampler in fieldnames(RLSSolver)
    @test :routine in fieldnames(RLSSolver)
    @test :log in fieldnames(RLSSolver)
    @test :stop in fieldnames(RLSSolver)
    @test :x in fieldnames(RLSSolver)

end

end # End module


text:
using StaticArrays, ArrayInterfaceCore, ArrayInterfaceStaticArraysCore, Test
using LinearAlgebra

x = @SVector [1,2,3]
@test ArrayInterfaceCore.ismutable(x) == false
@test ArrayInterfaceCore.ismutable(view(x, 1:2)) == false
@test ArrayInterfaceCore.can_setindex(typeof(x)) == false
@test ArrayInterfaceCore.buffer(x) == x.data

x = @MVector [1,2,3]
@test ArrayInterfaceCore.ismutable(x) == true
@test ArrayInterfaceCore.ismutable(view(x, 1:2)) == true

A = @SMatrix(randn(5, 5))
@test ArrayInterfaceCore.lu_instance(A) isa typeof(lu(A))
A = @MMatrix(randn(5, 5))
@test ArrayInterfaceCore.lu_instance(A) isa typeof(lu(A))

x = @SMatrix rand(Float32, 2, 2)
y = @SVector rand(4)
yr = ArrayInterfaceCore.restructure(x, y)
@test yr isa SMatrix{2, 2}
@test Base.size(yr) == (2,2)
@test vec(yr) == vec(Float32.(y))
z = rand(4)
zr = ArrayInterfaceCore.restructure(x, z)
@test zr isa SMatrix{2, 2}
@test Base.size(zr) == (2,2)
@test vec(zr) == vec(z)


text:
# automatically generated by the FlatBuffers compiler, do not modify

if !isdefined(@__MODULE__(), :MyGame) @__MODULE__().eval(:(module MyGame import FlatBuffers end)) end
if !isdefined(MyGame, :Sample) MyGame.eval(:(module Sample import FlatBuffers end)) end
include("MyGame/Sample/Color.jl")
include("MyGame/Sample/Weapon.jl")
include("MyGame/Sample/Equipment.jl")
include("MyGame/Sample/Vec3.jl")
include("MyGame/Sample/Monster.jl")

max_stars_repo_name: rjkat/flatbuffers-julia

max_issues_repo_name: rjkat/flatbuffers-julia

max_forks_repo_name: rjkat/flatbuffers-julia


text:
"""
# MomentTensors

The MomentTensors module provides routines for dealing with moment tensors,
such as creating them from double-couple parameters (strike, dip, rake, moment)
or evaluating source radiation patterns from them.

This module uses the Harvard/Global CMT convention for all input/output:
  - x // r         (local radial = upwards)
  - y // θ (theta) (local south = colatitude)
  - z // ϕ (phi)   (local east = (co)longitude)

(Individual routines may use different conventions internally.)

Other conventions are:

* Aki & Richards, Kennett:
  - x // N
  - y // E
  - z // down

MTs are stored internally as immutable types containing only a 6-element array
with the following construction:

    M = [Mrr, Mtt, Mpp, Mrt, Mrp, Mtp]
      = [Mrr, Mθθ, Mϕϕ, Mrθ, Mrϕ, Mθϕ]

They can then be accessed using vector index parameters or symbols like so:

    M[:tt] # = 1.5.
    mtp = M[:tp]
    # The following are all true
    M[2,3] == M[:θϕ]
    M[1,2] == Mrθ
    M[1,1] == M.m[1]

All angles are in degrees, always.
"""
module MomentTensors

using LinearAlgebra
using Rotations, StaticArrays

export
    MT,
    amplitude_v_azimuth,
    cmtsolution,
    decompose,
    eps_non_dc,
    m0,
    mw,
    ndk,
    radiation_pattern,
    rotate

"Array to convert from the 6-vector into a 3x3 tensor.  `m.m[_ij2k[i,j]] === m.m[k]`"
const _ij2k = @SMatrix [1 4 5
                        4 2 6
                        5 6 3]

"""
    MT(rr, θθ, ϕϕ, rθ, rϕ, θϕ) -> ::MT
    MT(M::AbstractVector(6)) -> ::MT
    MT(M::AbstractArray(3,3), warn=true) -> ::MT
    MT(strike, dip, rake, M0) -> ::MT

Construct a new `MT` (moment tensor) in the native frame used by `MomentTensors`:
- Radial (`r`) upwards
- Colatitude (`θ` or `t`) southwards
- Longitude (`ϕ` or `p`) eastwards)

Several forms exist to construct a moment tensor:
* Supply individual components as a list of arguments
* Supply a 6-vector
* Specify a strike, dip and rake in degrees, and scalar moment (N.m)
* Give a 3×3 matrix

For the last case, the upper half of the matrix only is used.  By default
symmetry is checked and a warning is issued if the input is not approximately
symmetric.  This can be turned off when `warn` is `false.

One may access the values of a moment tensor `M` in two ways:
1. M[i,j] yields the elements of `M.m` as if they were a two-tensor
2. M[::Symbol] yields the elements by name; see [`getindex`](@ref Base.getindex(::MT, ::Symbol)) for details
"""
struct MT{T<:Number}
    m::SVector{6,T}
    MT{T}(m) where T = new{T}(m)
end

MT(args...) = MT{float(promote_type(eltype.(args)...))}(args...)
MT(m::Union{AbstractVector{T}, AbstractArray{T,2}}) where T = MT{float(T)}(m)
MT{T}(rr, tt, pp, rt, rp, tp) where T = MT{T}(SVector{6,T}(rr, tt, pp, rt, rp, tp))
MT{T}(strike, dip, rake, M0) where T = MT{T}(_sdr2mt(strike, dip, rake, M0))
MT(strike, dip, rake, M0) = MT{float(typeof(M0))}(_sdr2mt(strike, dip, rake, M0))
function MT{T}(m::AbstractArray{U,2} where U, warn=true) where T
    size(m) == (3,3) ||
        throw(ArgumentError("2-dimensional array must have dimensions `(3,3)` for a `MT`"))
    warn && !(m[1,2] ≈ m[2,1] && m[1,3] ≈ m[3,1] && m[2,3] ≈ m[3,2]) &&
        @warn("Supplied tensor is not symmetric; using upper elements")
    MT{T}(SVector{6,T}(m[1,1], m[2,2], m[3,3], m[1,2], m[1,3], m[2,3]))
end

"""
    _matrix(m::MT)

Return a 3×3 `SMatrix` for the moment tensor `m`.
"""
_matrix(m::MT) = SMatrix{3,3}(m[1], m[4], m[5],
                              m[4], m[2], m[6],
                              m[5], m[6], m[3])

# Overloaded base operators and constructors
Base.Array(m::MT) = Array(_matrix(m))

"""
    getindex(m::MT, i, j) -> val
    getindex(m::MT, s::Symbol) -> val

Return a component of the moment tensor, `m`.

Two forms are permitted.  Either two indices may be supplied, `i` and `j`, or
a symbol giving the name of the component.  This must be one of:
    :rr, :tt, :pp, :rt, :tr, :rp, :pr, :tp, :pt
"""
Base.getindex(m::MT, s::Symbol) = m.m[_symbol2index(s)]
Base.getindex(m::MT, i) = m.m[i]
Base.getindex(m::MT, i, j) = m.m[_ij2k[i,j]]
Base.getindex(m::MT, i, j, inds...) = m[inds[i],inds[j]][inds...]
Base.size(m::MT) = (3, 3)
Base.isapprox(m1::MT, m2::MT; kwargs...) = isapprox(m1.m, m2.m; kwargs...)
Base.eltype(::MT{T}) where T = T

Base.:+(m::MT, a::Number) = MT(m.m .+ a)
Base.:+(a::Number, m::MT) = MT(a .+ m.m)
Base.:+(a::MT, b::MT) = MT(a.m .+ b.m)
Base.:-(m::MT, a::Number) = MT(m.m .- a)
Base.:-(a::Number, m::MT) = MT(a .- m.m)
Base.:-(a::MT, b::MT) = MT(a.m .- b.m)
Base.:-(m::MT) = MT(-m.m)
Base.:*(m::MT, a::Number) = MT(m.m.*a)
Base.:*(a::Number, m::MT) = MT(a.*m.m)
Base.:/(m::MT, a::Number) = MT(m.m./a)
Base.:/(a::Number, m::MT) = MT(a./m.m)

# Routines using or manipulating MTs
"""
    amplitude_v_azimuth(m::MT, inclination, azimuth_range=0:359) -> P, SV, SH, j

Return the P, SV and SH wave amplitudes, `P`, `SV`, `SH` respectively, and source
polarisation, `j`, across the range of azimuth `azimuth_range`.
"""
function amplitude_v_azimuth(m::MT{T}, inclination, azimuth_range=0:359) where T
    p = Array{T}(undef, length(azimuth_range))
    sv = similar(p)
    sh = similar(p)
    j = Array{typeof(m[1]/oneunit(m[1]))}(undef, length(azimuth_range))
    for (i, az) in enumerate(azimuth_range)
        p[i], sv[i], sh[i], j[i] = radiation_pattern(m, az, inclination)
    end
    p, sv, sh, j
end

"""
    auxplane(strike, dip, rake) -> strike, dip, rake

Return the strike, dip and rake of the auxiliary plane, given the `strike`,
`dip` and `rake` of the fault plane, all in degrees.
"""
function auxplane(strike, dip, rake)
    # Formula taken from Shearer, Introduction to Seismology
    s1 = deg2rad(strike)
    d1 = deg2rad(dip)
    r1 = deg2rad(rake)

    d2 = acos(sin(r1)*sin(d1))

    sr2 = cos(d1)/sin(d2)
    cr2 = -sin(d1)*cos(r1)/sin(d2)
    r2 = atan(sr2, cr2)

    s12 = cos(r1)/sin(d2)
    c12 = -1/(tan(d1)*tan(d2))
    s2 = s1 - atan(s12, c12)

    strike2 = rad2deg(s2)
    dip2 = rad2deg(d2)
    rake2 = rad2deg(r2)

    if dip2 > 90.0
       strike2 = strike2 + 180.0
       dip2 = 180.0 - dip2
       rake2 = 360.0 - rake2
    end
    rake2 = mod(rake2 + 180.0, 360.0) - 180.0 # In range -180 to 180
    strike2 = mod(strike2, 360.0) # In range 0 to 360
    strike2, dip2, rake2
end

"""
    cmtsolution(str::String) -> m

Return the moment tensor `m` described by the contents of `str`, which are in
the CMTSOLUTION format used by SPECFEM3D(_GLOBE).
"""
function cmtsolution(str::AbstractString)
    lines = _getlines(str)
    length(lines) == 13 || throw(ArgumentError("The supplied string " *
        "does not have 13 lines.  Got:\n" * str))
    # Allow for Fortran double precision scientific notation (1d-3)
    vals = [parse(Float64, replace(split(l)[2], r"[dD]"=>"e")) for l in lines[8:13]]
    MT(vals.*1e-7) # Convert from dyne.cm to N.m
end

"""
    ndk(str::String) -> m, err

Return the moment tensor `m` and its uncertainty `err` given by the contents of
`str`, which are in the ndk format used by the Global CMT project.
"""
function ndk(str::AbstractString)
    lines = _getlines(str)
    length(lines) == 5 || throw(ArgumentError("The supplied string does " *
        "not have 5 lines.  Got:\n" * str))
    tokens = split(lines[4])
    length(tokens) == 13 || throw(ArgumentError("Line 4 of input does not " *
        "have 13 fields.  Got:\n" * lines[4]))
    expo = tryparse(Int, tokens[1])
    expo === nothing && throw(ArgumentError("First token of line 4 " *
        "of input is not an integer.  Got: '$(tokens[1])'"))
    expo = expo - 7 # Convert from dyne.cm to N.m
    m = MT(parse.(Float64, tokens[2:2:end]).*10.0^expo)
    err = MT(parse.(Float64, tokens[3:2:end]).*10.0^expo)
    m, err
end

_getlines(s::String) = split(chomp(s), '\n')

"""
    radiation_pattern(m::MT, azimuth, inclination) -> P, SV, SH, j

Calculate the far-field radiation pattern from the moment tensor `m`, along `azimuth`
measured from local north towards east, and `inclination`, measured away from
downwards towards the radial direction.  Angles in degrees.

### Diagram

    ^ V (SV)
    |
    |     /       View looking along ray
    |_ j /
    | ` /
    | |/
    | /
    |/
    x----------> H (SH)

### Returns

* P, SV, SH : Amplitudes of ray in P, SV and SH directionso
* j         : Source polarisation measured from SV towards SH (i.e., like ϕ′ in splitting)
"""
function radiation_pattern(m::MT{T}, azimuth, inclination) where T
    #=
    Uses formula given in pp. 70 ff. of

        The Seismic Wavefield, Volume 1.  Kennett, B.L.N., Cambridge University Press.

    Kennett uses:

    * x // N
    * y // E
    * z // down

    and:

    - i is angle from z towards x-y plane (incidence angle away from down)
    - phi is angle from x towards y (i.e., azimuth clocwise from N)
    - V is direction upwards (radial) when looking along ray at source
    - H is direction to the right when looking along ray at source:
    =#
    Mxx =  m[:tt]
    Myy =  m[:pp]
    Mzz =  m[:rr]
    Mxy = -m[:tp]
    Mxz =  m[:rt]
    Myz = -m[:rp]

    # Convert to radians in same precision as `m`, removing units
    angle_type = typeof(one(T))
    a = angle_type(deg2rad(azimuth))
    i = angle_type(deg2rad(inclination))

    # Some shortcuts
    sini = sin(i)
    cosi = cos(i)
    sinphi = sin(a)
    cosphi = cos(a)
    sintwophi = sin(2a)
    costwophi = cos(2a)
    # Calculate radiation pattern for P, SV and SH
    P = (sini^2)*(Mxx*cosphi^2 + Mxy*sintwophi + Myy*sinphi^2 - Mzz) +
        2*sini*cosi*(Mxz*cosphi + Myz*sinphi) + Mzz
    SV = sini*cosi*(Mxx*cosphi^2 + Mxy*sintwophi + Myy*sinphi^2 - Mzz) +
         cos(2i)*(Mxz*cosphi + Myz*sinphi)
    SH = sini*((Myy - Mxx)*sinphi*cosphi + Mxy*costwophi) +
         cosi*(Myz*cosphi - Mxz*sinphi)
    # Source polarisation in ray frame, measured from upwards towards the right
    j = rad2deg(atan(SV, SH))
    P, SV, SH, j
end

"""
    rotate(m::MT, angle_r, angle_θ, angle_ϕ) -> ::MT

Return a rotated version of a moment tensor, rotated in turn about the axes
r, θ and ϕ.  The rotation appears clockwise when looking down each axis
towards the origin.
"""
function rotate(m::MT{T}, r, t, p) where T
    R = RotZ(deg2rad(p)) * RotY(deg2rad(t)) * RotX(deg2rad(r))
    m′ = R'*_matrix(m)*R
    MT{T}(m′)
end

"""
    mw(m0) -> mw

Convert a scalar moment (Nm) to a moment magnitude
"""
mw(M0) = 2/3*log10(M0*1e7/oneunit(M0)) - 10.7

"""
    m0(mw) -> m0

Convert a moment magnitude to a scalar moment (Nm)
"""
m0(Mw) = 10^(3/2*(Mw + 10.7))/1e7

"""
    _sdr2mt(strike, dip, rake, m0) -> ::SVector{6}

Return a moment tensor defined by the `strike`, `dip` and `rake` of the fault
plane, and a scalar moment `m0`.  Use `m0(mw)` if you want to specify the
moment mangitude instead.

Convention:
    * Strike is measured such that the downdip direction is 90° greater in azimuth
      than the strike.
    * Dip is measured between 0 and 90°.
    * Rake is measured from the strike direction, positive anticlockwise on the
      footwall, such that positive rakes imply a thrust component, and negative
      rakes a normal component.  Rakes should be in the range -180° to 180°

All angles are in degrees.
"""
function _sdr2mt(strike, dip, rake, M0)
    0 <= dip <= 90 || throw(ArgumentError("`dip` must be between 0° and 90°"))
    s, d, r = deg2rad(strike), deg2rad(dip), deg2rad(rake)
    tt = -M0*(sin(d)*cos(r)*sin(2s) + sin(2d)*sin(r)*sin(s)^2)
    tp = -M0*(sin(d)*cos(r)*cos(2s) + sin(2d)*sin(r)*sin(2s)/2)
    rt = -M0*(cos(d)*cos(r)*cos(s)  + cos(2d)*sin(r)*sin(s))
    pp =  M0*(sin(d)*cos(r)*sin(2s) - sin(2d)*sin(r)*cos(s)^2)
    rp =  M0*(cos(d)*cos(r)*sin(s)  - cos(2d)*sin(r)*cos(s))
    rr =  M0*sin(2d)*sin(r)
    @SVector [rr, tt, pp, rt, rp, tp]
end

"""
    _symbol2index(s::Symbol) -> k

Return the index in the range `1 <= k <= 6` with which to access the `m` field
of an `MT` type.

The following are permitted; order does not matter:
    :rr, :tt, :θθ, :pp, :ϕϕ, :rt, :rθ, :rp, :rϕ, :tp, :θϕ
"""
function _symbol2index(s::Symbol)
    if     s == :rr                  1
    elseif s in (:tt, :θθ)           2
    elseif s in (:pp, :ϕϕ)           3
    elseif s in (:rt, :tr, :rθ, :θr) 4
    elseif s in (:rp, :pr, :rϕ, :ϕr) 5
    elseif s in (:tp, :pt, :θϕ, :ϕθ) 6
    else throw(ArgumentError("'$s' is not a valid MT component name"))
    end
end

"""
    eps_non_dc(m::MT) -> ϵ

Return the deviatoric non-double-couple component `ϵ` of the moment tensor `m`,
defined by Giardini (1983) as

ϵ = –λ₂/max(|λ₁|, |λ₃|),

where the eigenvalues of the moment tensor are λ₁ ≥ λ₂ ≥ λ₃.

## References
- Giardini, D., 1983. Regional deviation of earthquake source mechanisms from the
  'double couple' model. In: H. Kanamori and E. Boschi (Editors), Earthquakes:
  Observation North-Holland, Amsterdam, Theory and Interpretation, pp. 345-353.
"""
function eps_non_dc(m::MT{T}) where T
    # Remove units if present
    m′ = m/oneunit(T)
    # eigen returns eigenvalues sorted by default
    λ = eigen(Symmetric(_matrix(m′))).values
    -λ[2]/max(abs(λ[1]), abs(λ[3]))
end

# Printing to the REPL looks like a 3×3 matrix
function Base.show(io::IO, mime::MIME"text/plain", m::MT{T}) where T
    println(io, "MomentTensors.MT{$T}:")
    io = IOContext(io, :compact=>true)
    Base.print_array(io, Array(m))
end

"""
    MTDecomposition{T<:Number}

Struct containing the results of [`decompose`](@ref).

Accessible fields are:

- `isotropic`: isotropic part of the moment tensor
- `isotropic_moment`: isotropic moment
- `isotropic_prop`: isotropic proportion of moment tensor
- `deviatoric`: deviatoric part of the moment tensor
- `deviatoric_moment`: deviatoric moment
- `deviatoric_prop`: deviatoric_prop part of the moment tensor
- `clvd`: CLVD part of the moment tensor
- `clvd_moment`: CLVD moment
- `clvd_prop`: CLVD proportion of moment tensor
- `double_couple`: double couple part of the moment tensor
- `double_couple_moment`: double couple moment
- `double_couple_prop`: double couple proportion of moment tensor
"""
struct MTDecomposition{T<:Number}
    isotropic::MT{T}
    isotropic_moment::T
    isotropic_prop::T
    deviatoric::MT{T}
    deviatoric_moment::T
    deviatoric_prop::T
    clvd::MT{T}
    clvd_moment::T
    clvd_prop::T
    double_couple::MT{T}
    double_couple_moment::T
    double_couple_prop::T
end

"""
    decompose(m::MT) -> iso, dev, dc, clvd, iso_m0, dev_m0, prop_iso, prop_dev, prop_dc, prop_clvd, m0

Decompose the arbitrary moment tensor `m` into its isotropic, `iso`,
and deviatoric parts, `dev`, plus the CLVD, `clvd`, and double-couple,
`dc`, components.  The proportion of the moment tensor represented by
the isotropic, deviatoric, double-couple and CLVD components are
given respectively by `prop_iso`, `prop_dev`, `prop_dc` and `prop_clvd`.
The scalar moment in Nm is given by `m0`.
"""
function decompose(m::MT{T}) where T
    # Remove units if present for eigendecomposition
    m = m/oneunit(T)
    @inbounds tr = m[1] + m[2] + m[3]
    _zero = zero(one(T))
    m_iso = MT(tr/3, tr/3, tr/3, _zero, _zero, _zero) # Isotropic MT
    m0_iso = abs(tr)/3 # Moment of isotropic part
    m_dev = m - m_iso # Deviatoric MT
    evals, evecs = LinearAlgebra.eigen(Symmetric(_matrix(m_dev)))
    # Sort eigenvectors (basis vectors) by absolute eigenvalues
    i = _sortperm_abs_3(evals)
    @inbounds begin
        e1 = SVector{3}(evecs[:,i[1]])
        e2 = SVector{3}(evecs[:,i[2]])
        e3 = SVector{3}(evecs[:,i[3]])
    end
    # λ₁ is smallest absolute eigenvalue, corresponding to e1
    @inbounds λ₁, λ₂, λ₃ = evals[i]
    # Largest absolute eigenvector gives moment of deviatoric part
    m0_dev = abs(λ₃)
    # Proportion of
    F = m0_dev < eps(T) ? T(0.5)/oneunit(T) : -λ₁/λ₃
    m_dc = MT(λ₃*(1 - 2F)*(e3.*e3' - e2.*e2'))
    m_clvd = m_dev - m_dc
    m0 = m0_iso + m0_dev
    iso_prop = m0_iso/m0
    dev_prop = m0_dev/m0
    dc_prop = (1 - 2F)*(1 - iso_prop)
    clvd_prop = 1 - iso_prop - dc_prop
    # Return units if any
    u = oneunit(T)
    return (iso=u*m_iso, dev=u*m_dev, dc=u*m_dc, clvd=u*m_clvd, iso_m0=u*m0_iso, dev_m0=u*m0_dev,
        prop_iso=iso_prop, prop_dev=dev_prop, prop_dc=dc_prop, prop_clvd=clvd_prop,
        m0=u*m0)
end

# Taken from StaticArrays.jl
#    https://github.com/JuliaArrays/StaticArrays.jl/blob/a0179213b741c0feebd2fc6a1101a7358a90caed/src/eigen.jl#L267
# which is distributed under the MIT "Expat" licence.
"""
    _sortperm_abs_3(vector) -> perm

Given a length-3 `vector`, return the permutation vector which sorts the
values in increasing absolute value.  `perm` is a length-3 `SVector`.

(Note that if `vector` is longer than 3, then elements after the first
three are ignored and `perm` simply gives the permutations of the
first three elements.)
"""
@inline function _sortperm_abs_3(v)
    abs_v = abs.(v)
    perm = SVector(1, 2, 3)
    (abs_v[perm[1]] > abs_v[perm[2]]) && (perm = SVector(perm[2], perm[1], perm[3]))
    (abs_v[perm[2]] > abs_v[perm[3]]) && (perm = SVector(perm[1], perm[3], perm[2]))
    (abs_v[perm[1]] > abs_v[perm[2]]) && (perm = SVector(perm[2], perm[1], perm[3]))
    perm
end

end # module


text:
# Part One: using direct collection method
using InfiniteOpt, Ipopt, Distributions, Plots

# Set the SEIR parameters
b = 0.525
d = 0.5
c = 0.0001
e = 0.5
g = 0.1
a = 0.2
A = 0.1

# Set the domain information
t0 = 0
tf = 20

# Set the intial condition values
s0 = 1000
e0 = 100
i0 = 50
n0 = 1165

model = InfiniteModel(Ipopt.Optimizer)

@infinite_parameter(model, t ∈ [t0, tf], num_supports = 101,
    derivative_method = OrthogonalCollocation(2))

@variable(model, S ≥ 0, Infinite(t))
@variable(model, E ≥ 0, Infinite(t))
@variable(model, I ≥ 0, Infinite(t))
@variable(model, N ≥ 0, Infinite(t))
@variable(model, 0 ≤ u ≤ 0.9, Infinite(t), start = 0.2)

@objective(model, Min, ∫(A * I + u^2, t))

# Define the initial conditions
@constraint(model, S(0) == s0)
@constraint(model, E(0) == e0)
@constraint(model, I(0) == i0)
@constraint(model, N(0) == n0)

# Define the SEIR equations
@constraint(model, s_constr, ∂(S, t) == b * N - d * S - c * S * I - u * S)
@constraint(model, e_constr, ∂(E, t) == c * S * I - (e + d) * E)
@constraint(model, i_constr, ∂(I, t) == e * E - (g + a + d) * I)
@constraint(model, n_constr, ∂(N, t) == (b - d) * N - a * I)

# Define the infection rate limit

print(model)

optimize!(model)

I_value = value(I)
u_value = value(u)
tspan = value(t)
plot(tspan, I_value)
plot(tspan, u_value)

# Part Two: using Deep Learning Method
using DiffEqFlux, DifferentialEquations
const solveeq = DifferentialEquations.solve

tspan = (0.0f0, 20.0f0)
ann = FastChain(FastDense(1, 32, tanh), FastDense(32, 1), (x, p) -> 0.9 * sigmoid.(x))
θ = initial_params(ann)
function sir_nn(du, u, p, t)
    S, E, I, N = u
    du[1] = b * N - d * S - c * S * I - ann([t], p)[1] * S
    du[2] = c * S * I - (e + d) * E
    du[3] = e * E - (g + a + d) * I
    du[4] = (b - d) * N - a * I
end
u0 = [s0, e0, i0, n0]
ts = Float32.(collect(0.0:0.2:tspan[2]))
prob = ODEProblem(sir_nn, u0, tspan, θ)
sol_init = solveeq(prob, Vern9(), abstol = 1e-10, reltol = 1e-10)
plot(sol_init)
function predict_adjoint(θ)
    Array(solveeq(prob, Vern9(), p = θ, saveat = ts, sensealg = InterpolatingAdjoint(autojacvec = ReverseDiffVJP(true))))
end
function loss_adjoint(θ)
    x = predict_adjoint(θ)
    mean(abs, A .* x[3, :]) + mean(abs2, [first(ann([t], θ)) for t in ts])
end
l = loss_adjoint(θ)
cb = function (θ, l)
    println(l)
    p = plot(solveeq(remake(prob, p = θ), Tsit5(), saveat = 0.2), lw = 3)
    #plot!(p, ts, [first(ann([t], θ)) for t in ts], label = "u(t)", lw = 3)
    display(p)
    return false
end
# Display the ODE with the current parameter values.
cb(θ, l)
loss1 = loss_adjoint(θ)
res1 = DiffEqFlux.sciml_train(loss_adjoint, θ, ADAM(0.05), cb = cb, maxiters = 100)
res2 = DiffEqFlux.sciml_train(loss_adjoint, res1.u,
    BFGS(initial_stepnorm = 0.01), cb = cb, maxiters = 100,
    allow_f_increases = false)

# Part Three: Saving the results
using DataFrames
using CSV
results_saving_opt = DataFrame()
results_saving_opt.t = tspan
results_saving_opt.I = I_value
results_saving_opt.u = u_value
results_saving_opt.I_nn = Array(solveeq(remake(prob, p = res2.u), Tsit5(), saveat = 0.2))[3, :]
results_saving_opt.u_nn = [first(ann([t], re2.u)) for t in ts]
CSV.write("results_saving_SEINmodel.csv", results_saving_opt)

plot(results_saving_opt.t, results_saving_opt.I, lw = 2, label = "I(t) by direct collection")
plot!(results_saving_opt.t, results_saving_opt.I_nn, lw = 2, label = "I(t) by deep learning")
xlabel!("t(day)")
ylabel!("I(t)")
savefig("sein_I.png")

plot(results_saving_opt.t[2:end-1], results_saving_opt.u[2:end-1], lw = 2, label = "u(t) by direct collection")
plot!(results_saving_opt.t[2:end-1], results_saving_opt.u_nn[2:end-1], lw = 2, label = "u(t) by deep learning")
xlabel!("t(day)")
ylabel!("u(t)")
savefig("sein_control.png")

# Loss saving

opt_method_loss = 13.90369069230782
nn_method_loss = 13.711496454176164