text:
using Pumas, PumasTutorials, CSV, Plots

data = PumasTutorials.tutorial_data("data/nca","SD_IV_infusion")
data = CSV.read(data,missingstring="NA")
first(data,10)

timeu = u"hr"
concu = u"mg/L"
amtu  = u"mg"

pop = read_nca(data, id=:ID, time=:time, conc=:DV, amt=:DOSE, ii=24timeu,
    route=:Formulation, duration=:Infusion_Time,timeu=timeu,
    concu=concu, amtu=amtu,lloq=0.4concu)

plot(pop)

NCA.auc(pop,auctype=:last,method=:linear)

NCA.auc(pop,auctype=:inf,method=:linuplogdown)

NCA.auc(pop, interval=(0,12).*timeu)

NCA.auc(pop, interval=[(0,12).*timeu,(0,6).*timeu])

NCA.lambdaz(pop)

NCA.lambdaz(pop, threshold=3)
NCA.lambdaz(pop, idxs=[18,19,20])
NCA.lambdaz(pop, slopetimes=[18.5,19,19.5].*timeu)

cmax = NCA.cmax(pop[1])

 NCA.normalizedose(cmax,pop[1])

NCA.lambdazr2(pop)
NCA.lambdazadjr2(pop)
NCA.lambdazintercept(pop)
NCA.lambdaztimefirst(pop)
NCA.lambdaznpoints(pop)

NCA.tmax(pop)
NCA.cmin(pop)
NCA.tmin(pop)

NCA.tlast(pop)
NCA.clast(pop)

NCA.aumc(pop)
NCA.aumclast(pop)

NCA.thalf(pop)

NCA.cl(pop)

NCA.vss(pop)
NCA.vz(pop)

report = NCAReport(pop)
report = NCA.to_dataframe(report)

report = NCAReport(pop,pred=true)
report = NCA.to_dataframe(report)

CSV.write("./tutorials/nca/report_SD_IV_infusion.csv", report)

using PumasTutorials
PumasTutorials.tutorial_footer(WEAVE_ARGS[:folder],WEAVE_ARGS[:file])


text:
# License for this file: MIT (expat)
# Copyright 2017-2018, DLR Institute of System Dynamics and Control
#
# This file is part of module
#   Modia3D.Shapes (Modia3D/Solids/_module.jl)
#

# MassProperties consisting of mass m, center of mass rCM and inertia tensor I
struct MassProperties <: Modia3D.AbstractMassProperties
    m::Float64                 # mass [kg]
    rCM::Frames.Vector3D       # position vector from object3d frame to center of mass resolved in object3d frame [m]
    I::SMatrix{3,3,Float64,9}  # inertia matrix w.r.t. center of mass resolved in object3d frame [kg.m^2]

    #---------------- different constructors for MassProperties -----------------
    # Constructor 0: takes mass, centerOfMass and inertiaMatrix as input values
    function MassProperties(mass::Number, centerOfMass::AbstractVector, inertiaMatrix::AbstractMatrix)
        @assert(mass > 0.0)
        new(mass, centerOfMass, inertiaMatrix)
    end
end

struct MassPropertiesFromShape <: Modia3D.AbstractMassProperties
    function MassPropertiesFromShape()
        new()
    end
end

struct MassPropertiesFromShapeAndMass <: Modia3D.AbstractMassProperties
    mass::Number   # mass in [kg]

    function MassPropertiesFromShapeAndMass(;mass::Number=1.0)
        @assert(mass >= 0.0)
        new(mass)
    end
end

function Base.show(io::IO, mp::MassProperties)
    print(io, "mass = ", mp.m,
            ", centerOfMass = ", mp.rCM,
            ", Ixx = ", mp.I[1,1],
            ", Iyy = ", mp.I[2,2],
            ", Izz = ", mp.I[3,3],
            ", Ixy = ", mp.I[1,2],
            ", Ixz = ", mp.I[1,3],
            ", Iyz = ", mp.I[2,3])
end

# Constructor a: mass, centerOfMass and entries of inertia tensor are optional
#                --> if nothing special is defined it takes predefined values (= zero values)
MassProperties(; mass::Number=0.0, centerOfMass=Modia3D.ZeroVector3D,
               Ixx::Number=0.0, Iyy::Number=0.0, Izz::Number=0.0,
               Ixy::Number=0.0, Ixz::Number=0.0, Iyz::Number=0.0) =
                  MassProperties(mass, centerOfMass, [Ixx Ixy Ixz; Ixy Iyy Iyz; Ixz Iyz Izz])
# Constructor b: shape and mass is given, center of mass and inertia tensor is
#                calculated via shape --> constructor 0 is called
MassProperties(shape::Modia3D.AbstractGeometry, mass::Number) =
                     MassProperties(mass, centroid(shape), inertiaMatrix(shape,mass))
# Constructor c: shape and material is given, mass is computed via volume of
#                shape and density --> constructor b is called
MassProperties(shape::Modia3D.AbstractGeometry, material::SolidMaterial) =
                     MassProperties(shape, material.density*volume(shape))
# Constructor d: shape and materialName is given, material must be defined in
#                solidMaterialPalette --> constructor c is called
MassProperties(shape::Modia3D.AbstractGeometry, materialName::AbstractString) =
                     MassProperties(shape, solidMaterialPalette[materialName])

# structure InternalMassProperties is only for internal purposes
mutable struct InternalMassProperties <: Modia3D.AbstractMassProperties
   m::Float64                 # mass in [kg]
   rCM::Frames.Vector3D   # center of mass in [m]
   I::SMatrix{3,3,Float64,9}  # inertia matrix in [kg.m^2]

   InternalMassProperties() = new(0.0, Modia3D.ZeroVector3D,
                                  SMatrix{3,3,Float64,9}(0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0))
end

#=
getMassProperties(massProperties::Union{Modia3D.AbstractMassProperties, Number, AbstractString, SolidMaterial, Nothing},
                  shape::Modia3D.AbstractGeometry, solidMaterial::Union{AbstractString,SolidMaterial,Nothing}) = createMassProperties(massProperties, shape, solidMaterial)
=#
createMassProperties(massProperties::MassProperties, shape, solidMaterial) = massProperties

createMassProperties(massProperties::Union{Number, SolidMaterial, AbstractString}, shape::Modia3D.AbstractGeometry, solidMaterial) = MassProperties(shape, massProperties)

# compute mass properties from shape and material
function createMassProperties(massProperties::Union{MassPropertiesFromShape,Nothing}, shape::Modia3D.AbstractGeometry, solidMaterial::Union{AbstractString,SolidMaterial,Nothing})
    if isnothing(solidMaterial)
        error("It is not possible to compute mass properties (MassPropertiesFromShape = ", massProperties,") for shape = ", shape , " because no solidMaterial is defined.")
    else
        return Modia3D.MassProperties(shape, solidMaterial)
    end
end

# compute mass properties from shape and mass
function createMassProperties(massProperties::MassPropertiesFromShapeAndMass, shape::Modia3D.AbstractGeometry, solidMaterial::Union{AbstractString,SolidMaterial,Nothing})
    return Modia3D.MassProperties(shape, massProperties.mass)
end

function createMassProperties(massProperties::MassPropertiesFromShape, shape::Nothing, solidMaterial::Union{AbstractString,SolidMaterial,Nothing})
    error("It is not possible to compute mass properties (MassPropertiesFromShape = ", massProperties,") because no shape is defined.")
end

function createMassProperties(massProperties::MassPropertiesFromShapeAndMass, shape::Nothing, solidMaterial::Union{AbstractString,SolidMaterial,Nothing})
    error("It is not possible to compute mass properties (MassPropertiesFromShapeAndMass) if no shape is defined.")
end

ext: jl


text:
using Genser
using Test
import JSON

@testset "Genser.jl" begin
    include("datamodel_tests.jl")
    include("derivetype_tests.jl")
    include("convert_tests.jl")
    include("registry_tests.jl")
    include("serialize_tests.jl")
    include("deserialize_tests.jl")
end


text:
using Test

using JuMP
using MIPVerify: owner_model, interval, lower_bound, upper_bound
@isdefined(TestHelpers) || include("../TestHelpers.jl")

TestHelpers.@timed_testset "ConditionalJuMP.jl" begin
    @testset "owner_model" begin
        m = Model()
        x = @variable(m)
        @test owner_model(x) == m

        y = @variable(m)
        @test owner_model(3 * x + 5 * y) == m
        @test owner_model([x, y]) == m
        @test owner_model([2 * x + 3 * y, 3 * x + 5 * y]) == m
    end

    @testset "interval, lower_bound, upper_bound" begin
        m = Model()
        @variable(m, 1 <= x <= 3)
        @variable(m, 2 <= y <= 5)
        e = 2 * x + 1

        @testset "Number" begin
            i = interval(5.0)
            @test lower_bound(i) == 5.0
            @test upper_bound(i) == 5.0
        end

        @testset "JuMP.VariableRef" begin
            @test lower_bound(x) == 1
            @test upper_bound(x) == 3

            i = interval(x)
            @test lower_bound(i) == 1
            @test upper_bound(i) == 3
        end

        @testset "JuMP.GenericAffExpr" begin
            @test lower_bound(e) == 3
            @test upper_bound(e) == 7

            @test lower_bound(e - 2x) == 1
            @test upper_bound(e - 2x) == 1

            @test lower_bound(2x - 3y) == -13
            @test upper_bound(2x - 3y) == 0

            @test lower_bound(AffExpr(2)) == 2
            @test upper_bound(AffExpr(2)) == 2

            i = interval(e)
            @test lower_bound(i) == lower_bound(e)
            @test upper_bound(i) == upper_bound(e)
        end

        @testset "Interval" begin
            i = interval(5.0)
            @test lower_bound(i) == 5.0
            @test upper_bound(i) == 5.0
        end
    end
end


text:
include("../src/TSne.jl")
using RDatasets
using Gadfly
using TSne

lables = ()

println("Using Iris dataset.")
iris = dataset("datasets", "iris")
X = array(iris[:,1:4])
labels = iris[:,5]
plotname = "iris"
initial_dims = -1
iterations = 1500
perplexity = 15

println("X dimensions are: " * string(size(X)))
Y = tsne(X, 2, initial_dims, iterations, perplexity)
println("Y dimensions are: " * string(size(Y)))

theplot = plot(x=Y[:,1], y=Y[:,2], color=labels)

draw(PDF(plotname*".pdf", 4inch, 3inch), theplot)
draw(SVG(plotname*".svg", 4inch, 3inch), theplot)


text:
istransparent(::Any) = false
istransparent(::MLJType) = true

"""
    params(m::MLJType)

Recursively convert any transparent object `m` into a named tuple,
keyed on the fields of `m`. An object is *transparent* if
`MLJModelInterface.istransparent(m) == true`. The named tuple is
possibly nested because `params` is recursively applied to the field
values, which themselves might be transparent.

Most objects of type `MLJType` are transparent.

```julia
julia> params(EnsembleModel(atom=ConstantClassifier()))
(atom = (target_type = Bool,),
weights = Float64[],
bagging_fraction = 0.8,
rng_seed = 0,
n = 100,
parallel = true,)
```
"""
params(m) = params(m, Val(istransparent(m)))
params(m, ::Val{false}) = m

function params(m, ::Val{true})
    fields = fieldnames(typeof(m))
    return NamedTuple{fields}(Tuple([params(getfield(m, field)) for field in fields]))
end


text:
struct GIValue{G <: AbstractGrid}
    grid::G
    gdata::Vector{Float64}
end

evaluate(v::GIValue, s::AbstractVector{Float64}) = interpolate(v.grid, v.gdata, convert(Vector{Float64}, s))

@with_kw mutable struct CWorldSolver{G<:AbstractGrid, RNG<:AbstractRNG} <: Solver
    grid::G                     = RectangleGrid(range(0.0, stop=10.0, length=30), range(0.0, stop=10.0, length=30))
    max_iters::Int              = 50
    tol::Float64                = 0.01
    m::Int                      = 20
    value_hist::AbstractVector  = []
    rng::RNG                    = Random.GLOBAL_RNG
end

struct CWorldPolicy{V} <: Policy
    actions::Vector{Vec2}
    Qs::Vector{V}
end

function POMDPs.solve(sol::CWorldSolver, w::CWorld)
    sol.value_hist = []
    data = zeros(length(sol.grid))
    val = GIValue(sol.grid, data)

    for k in 1:sol.max_iters
        newdata = similar(data)
        for i in 1:length(sol.grid)
            s = Vec2(ind2x(sol.grid, i))
            if isterminal(w, s)
                newdata[i] = 0.0
            else
                best_Qsum = -Inf
                for a in actions(w, s)
                    Qsum = 0.0
                    for j in 1:sol.m
                        sp, r = @gen(:sp, :r)(w, s, a, sol.rng)
                        Qsum += r + discount(w)*evaluate(val, sp)
                    end
                    best_Qsum = max(best_Qsum, Qsum)
                end
                newdata[i] = best_Qsum/sol.m
            end
        end
        push!(sol.value_hist, val)
        print("\rfinished iteration $k")
        val = GIValue(sol.grid, newdata)
    end

    print("\nextracting policy...     ")

    acts = collect(actions(w))
    Qs = Vector{GIValue}(undef,length(acts))
    for j in 1:length(acts)
        a = acts[j]
        qdata = similar(val.gdata)
        for i in 1:length(sol.grid)
            s = Vec2(ind2x(sol.grid, i))
            if isterminal(w, s)
                qdata[i] = 0.0
            else
                Qsum = 0.0
                for k in 1:sol.m
                    sp, r = @gen(:sp, :r)(w, s, a, sol.rng)
                    Qsum += r + discount(w)*evaluate(val, sp)
                end
                qdata[i] = Qsum/sol.m
            end
        end
        Qs[j] = GIValue(sol.grid, qdata)
    end
    println("done.")

    return CWorldPolicy(acts, Qs)
end

function POMDPs.action(p::CWorldPolicy, s::AbstractVector{Float64})
    best = action_ind(p, s)
    return p.actions[best]
end

action_ind(p::CWorldPolicy, s::AbstractVector{Float64}) = argmax([evaluate(Q, s) for Q in p.Qs])

POMDPs.value(p::CWorldPolicy, s::AbstractVector{Float64}) = maximum([evaluate(Q, s) for Q in p.Qs])


text:
using Boh, Test, Luxor

mktempdir() do tmpdir
    cd(tmpdir) do
        Drawing(200, 200, "boh.svg")
        origin()
        boh()
        @test finish() == true
    end
end


text:
function Base.getindex(dbm::DBM, key::Union{String,Array{UInt8,1}})
  datum = gdbm_fetch(dbm.handle, Datum(key))
  value = unsafe_string(datum.dptr, datum.dsize)
  libc_free(datum.dptr)
  value
end

function Base.setindex!(dbm::DBM, value::Union{String,Array{UInt8,1}}, key::Union{String,Array{UInt8,1}})
  gdbm_store(dbm.handle, Datum(key), Datum(value))
  key
end

Base.start(dbm::DBM) = gdbm_firstkey(dbm.handle)

Base.done(dbm::DBM, state::Datum) = state.dptr == C_NULL

function Base.next(dbm::DBM, prev::Datum)
  key = unsafe_string(prev.dptr, prev.dsize)
  datum = gdbm_fetch(dbm.handle, prev)
  value = unsafe_string(datum.dptr, datum.dsize)
  libc_free(datum.dptr)
  next = gdbm_nextkey(dbm.handle, prev)
  libc_free(prev.dptr)
  ((key, value), next)
end

Base.length(dbm::DBM) = gdbm_count(dbm.handle)

Base.eltype(::Type{DBM}) = Union{Pair{String,String},Pair{String,Array{UInt8,1}},Pair{Array{UInt8,1},String},Pair{Array{UInt8,1},Array{UInt8,1}}}

Base.keytype(::Type{DBM}) = String

Base.valtype(::Type{DBM}) = String

Base.haskey(dbm::DBM, key::Union{String,Array{UInt8,1}}) = gdbm_exists(dbm.handle, Datum(key))

function Base.delete!(dbm::DBM, key::String)
  gdbm_delete(dbm.handle, Datum(key))
  dbm
end

function Base.pop!(dbm::DBM, key::Union{String,Array{UInt8,1}}, default = nothing)
  datum = Datum(key)
  if gdbm_exists(dbm.handle, datum)
    vdatum = gdbm_fetch(dbm.handle, datum)
    value = unsafe_string(vdatum.dptr, vdatum.dsize)
    libc_free(vdatum.dptr)
    gdbm_delete(dbm.handle, datum)
    return value
  end
  default isa String && return default
  throw(KeyError(key))
end

function Base.getkey(dbm::DBM, key::Union{String,Array{UInt8,1}}, default::String)
  gdbm_exists(dbm.handle, Datum(key)) && return key
  default
end

function Base.get(dbm::DBM, key::Union{String,Array{UInt8,1}}, default::String)
  value = default
  datum = Datum(key)
  if gdbm_exists(dbm.handle, datum)
    datum = gdbm_fetch(dbm.handle, datum)
    value = unsafe_string(datum.dptr, datum.dsize)
    libc_free(datum.dptr)
  end
  value
end

function Base.get!(dbm::DBM, key::Union{String,Array{UInt8,1}}, default::String)
  value = default
  datum = Datum(key)
  if gdbm_exists(dbm.handle, datum)
    vdatum = gdbm_fetch(dbm.handle, datum)
    value = unsafe_string(vdatum.dptr, vdatum.dsize)
    libc_free(vdatum.dptr)
  else
    gdbm_store(dbm.handle, datum, Datum(default))
  end
  value
end

function Base.isempty(dbm::DBM)
  gdbm_count(dbm.handle) == 0 && return true
  false
end

function Base.empty!(dbm::DBM)
  prev = gdbm_firstkey(dbm.handle)
  while prev.dptr ≠ C_NULL
    key = unsafe_string(prev.dptr, prev.dsize)
    next = gdbm_nextkey(dbm.handle, prev)
    gdbm_delete(dbm.handle, prev)
    libc_free(prev.dptr)
    prev = next
  end
end

function Base.in(pair::Union{Pair{String,String},Pair{String,Array{UInt8,1}},Pair{Array{UInt8,1},String},Pair{Array{UInt8,1},Array{UInt8,1}}}, dbm::DBM)
  datum = Datum(pair.first)
  !gdbm_exists(dbm.handle, datum) && return false
  vdatum = gdbm_fetch(dbm.handle, datum)
  value = unsafe_string(vdatum.dptr, vdatum.dsize)
  libc_free(vdatum.dptr)
  value == pair.second
end


text:
function PivotJoint(a::Body, b::Body, pivot::Vect)
    PivotJoint(ccall(dlsym(libchipmunk, :cpPivotJointNew), Ptr{Void}, (Ptr{Void}, Ptr{Void}, Vect,), a.ptr, b.ptr, pivot))
end

function PivotJoint(a::Body, b::Body, anchorA::Vect, anchorB::Vect)
    PivotJoint(ccall(dlsym(libchipmunk, :cpPivotJointNew2), Ptr{Void}, (Ptr{Void}, Ptr{Void}, Vect, Vect,), a.ptr, b.ptr, anchorA, anchorB))
end

function get_anchorA(joint::PivotJoint)
    ccall(dlsym(libchipmunk, :cpPivotJointGetAnchorA), Vect, (Ptr{Void},), joint.ptr)
end

function set_anchorA(joint::PivotJoint, anchorA::Vect)
    ccall(dlsym(libchipmunk, :cpPivotJointSetAnchorA), Void, (Ptr{Void}, Vect,), joint.ptr, anchorA)
end

function get_anchorB(joint::PivotJoint)
    ccall(dlsym(libchipmunk, :cpPivotJointGetAnchorB), Vect, (Ptr{Void},), joint.ptr)
end

function set_anchorB(joint::PivotJoint, anchorB::Vect)
    ccall(dlsym(libchipmunk, :cpPivotJointSetAnchorB), Void, (Ptr{Void}, Vect,), joint.ptr, anchorB)
end

export PivotJoint, get_anchorA, get_anchorB, set_anchorA, set_anchorB


text:
using ServerMetrics
using Test: @test, @test_throws
using Dates

# For testing purposes we also have MockStatsdBackend that simply stores the statds
# messages in a Set (unordered) and offers empty! method to clear the buffer.
Base.@kwdef mutable struct MockStatsdBackend <: ServerMetrics.AbstractServiceBackend
    messages::Set{String} = Set{String}()
end
ServerMetrics.send(b::MockStatsdBackend, msg::String) = push!(b.messages, msg)
function clear_messages!(b::MockStatsdBackend)
    b.messages = Set{String}()
end

# Every time we want to trigger statsd updates we want to clear existing
# messages from the MockStatsdBackend and we also want to wait for 1ms to
# ensure that the clock has time to advance.
# We have experienced issues where clock was running too slow and emissions
# and changes to metrics were happening at the same time.
function test_friendly_export_to_statsd(statsd::MockStatsdBackend, exporter::StatsdExporter)
    clear_messages!(statsd)
    # Ensure that the time of an emission is distinct from any other times
    # when metrics can be set by padding this call with small delays both
    # before and after.
    sleep(0.01)
    ServerMetrics.send_metric_updates(exporter)
    sleep(0.01)
end

@testcase "Counter changes" begin
    mock_statsd = MockStatsdBackend()
    r = MetricRegistry()
    c = register!(r, Counter(), "counter")
    exporter = StatsdExporter(metric_registries=Set([r]), statsd_backend=mock_statsd)

    test_friendly_export_to_statsd(mock_statsd, exporter)
    @test mock_statsd.messages == Set(["counter:0.0|c"])

    inc!(c, 1.0)
    test_friendly_export_to_statsd(mock_statsd, exporter)
    @test mock_statsd.messages == Set(["counter:1.0|c"])
    @test c.content.value[] == 1.0

    test_friendly_export_to_statsd(mock_statsd, exporter)
    @test mock_statsd.messages == Set()
    @test c.content.value[] == 1.0

    inc!(c, 2.0)
    test_friendly_export_to_statsd(mock_statsd, exporter)
    @test mock_statsd.messages == Set(["counter:2.0|c"])
    @test c.content.value[] == 3.0
end

@testcase "Gauge changes" begin
    mock_statsd = MockStatsdBackend()
    r = MetricRegistry()
    g = register!(r, Gauge(1.0), "gg")
    exporter = StatsdExporter(metric_registries=Set([r]), statsd_backend=mock_statsd)

    test_friendly_export_to_statsd(mock_statsd, exporter)
    @test mock_statsd.messages == Set(["gg:1.0|g"])

    test_friendly_export_to_statsd(mock_statsd, exporter)
    @test mock_statsd.messages == Set()

    inc!(g, 2.0)
    test_friendly_export_to_statsd(mock_statsd, exporter)
    @test mock_statsd.messages == Set(["gg:3.0|g"])

    dec!(g, 0.5)
    test_friendly_export_to_statsd(mock_statsd, exporter)
    @test mock_statsd.messages == Set(["gg:2.5|g"])

    inc!(g, 55.0)
    set!(g, 0.1)
    test_friendly_export_to_statsd(mock_statsd, exporter)
    @test mock_statsd.messages == Set(["gg:0.1|g"])
end

@testcase "Changes with labels" begin
    mock_statsd = MockStatsdBackend()
    r = MetricRegistry()
    reqs = register!(r, Counter(;action=String, response_code=Int64), "requests")
    inc!(reqs; action="get", response_code=200)
    inc!(reqs; action="put", response_code=500)
    inc!(reqs; action="put", response_code=500)

    exporter = StatsdExporter(metric_registries=Set([r]), statsd_backend=mock_statsd)
    test_friendly_export_to_statsd(mock_statsd, exporter)
    @test mock_statsd.messages == Set([
        "requests:1.0|c|#action:get,response_code:200",
        "requests:2.0|c|#action:put,response_code:500",
    ])
end

@testcase "Two metrics exported" begin
    mock_statsd = MockStatsdBackend()
    r = MetricRegistry()
    c = register!(r, Counter(), "cnt")
    g = register!(r, Gauge(1.0), "gg")
    exporter = StatsdExporter(metric_registries=Set([r]), statsd_backend=mock_statsd)

    test_friendly_export_to_statsd(mock_statsd, exporter)
    @test mock_statsd.messages == Set(["gg:1.0|g", "cnt:0.0|c"])

    inc!(c, 1.0)
    test_friendly_export_to_statsd(mock_statsd, exporter)
    @test mock_statsd.messages == Set(["cnt:1.0|c"])

    set!(g, 2.0)
    test_friendly_export_to_statsd(mock_statsd, exporter)
    @test mock_statsd.messages == Set(["gg:2.0|g"])

    inc!(c, 1.0)
    set!(g, 2.1)
    test_friendly_export_to_statsd(mock_statsd, exporter)
    @test mock_statsd.messages == Set(["cnt:1.0|c", "gg:2.1|g"])
end

@testcase "StatsdExporter with default registry" begin
    clear_registry!(get_default_registry())
    mock_statsd = MockStatsdBackend()
    r = MetricRegistry()
    c1 = register!(get_default_registry(), Counter(), "cnt_default_registry")
    c2 = register!(r, Counter(), "cnt_custom_registry")
    exporter = StatsdExporter(statsd_backend=mock_statsd)
    inc!(c1, 1.0)
    inc!(c2, 1.0)
    test_friendly_export_to_statsd(mock_statsd, exporter)
    @test mock_statsd.messages == Set(["cnt_default_registry:1.0|c"])
end

@testcase "StatsdExporter spawns background thread" begin
    mock_statsd = MockStatsdBackend()
    r = MetricRegistry()
    c = register!(r, Counter(), "cnt")
    exporter = StatsdExporter(
        send_interval=Millisecond(100),
        statsd_backend=mock_statsd,
        metric_registries=Set([r]))
    @test mock_statsd.messages == Set()
    inc!(c, 1.0)
    sleep(0.2)
    @test mock_statsd.messages == Set()
    start_statsd_exporter!(exporter)
    sleep(0.5)

    @test istaskstarted(exporter.periodic_task)
    @test !istaskdone(exporter.periodic_task)
    @test mock_statsd.messages == Set(["cnt:1.0|c"])
    t = stop_statsd_exporter!(exporter)
    @test isnothing(exporter.periodic_task)
    @test istaskdone(t)
end

@testcase "StatsdExporter not exporting with send_interval=0" begin
    mock_statsd = MockStatsdBackend()
    exporter = StatsdExporter(
        send_interval=Millisecond(0),
        statsd_backend=mock_statsd)
    start_statsd_exporter!(exporter)
    @test isnothing(exporter.periodic_task)
    @test isnothing(stop_statsd_exporter!(exporter))
end

# TODO(janrous): test that the emission generally does happen with
# roughly the expected frequency by waiting and counting number of
# emissions. reset! method may set the counters back to 0 for testing
# purposes.

@testcase "StatsdExporter sends UDP packets" begin
    # Create UDP listener and then send packet to it using send(::StatsdBackend, msg)
    server = UDPSocket()
    bind(server, ip"127.0.0.1", 12346)

    client = ServerMetrics.UDPBackend(ip"127.0.0.1", 12346)
    ServerMetrics.send(client, "test message")
    @test String(recv(server)) == "test message"
    close(server)
end


text:
include("defs.jl")

w = 0.2f0
bd = billiard_mushroom(1f0, w, 1f0, 0f0)

pc = MushroomTools.randomchaotic(1f0, w, 1f0)
pc = Particle(-0.01f0, 0.2f0, sqrt(3f0)/2)

pr = Particle(0.0f0, 1.2f0, 0.0f0)
pr2 = Particle(0.0f0, 1.9f0, 0.0f0)

ps = [pc, pr, pr2]
colors = [to_color(COLORS[i+1]) for i in 1:length(ps)]

fig, = interactive_billiard(bd, ps; colors = colors, tail = 10000)

# record_interaction(fig, joinpath(@__DIR__, "mushroom.mp4"))

max_stars_repo_name: JuliaDynamics/NonlinearDynamicsTextbook

max_issues_repo_name: JuliaDynamics/NonlinearDynamicsTextbook

max_forks_repo_name: JuliaDynamics/NonlinearDynamicsTextbook


text:
input_data = open("input03.txt") do file
    data = readlines(file)
    entries = map(x -> [String(y) for y in split(x, ",")], data)
    #numbers = map(x->parse(Int64, x), entries)
    #numbers
    entries
end

struct Point
    x::Int
    y::Int
end

function parse_prefix(path::Array{Point, 1}, inst::String) :: Array{Point, 1}
    prefix = inst[1]
    current_ending = path[end]
    distance = parse(Int64, inst[2:end])
    if prefix == 'R'
        push!(path, Point(current_ending.x + distance, current_ending.y))
    elseif prefix == 'L'
        push!(path, Point(current_ending.x - distance, current_ending.y))

    elseif prefix == 'D'
        push!(path, Point(current_ending.x, current_ending.y + distance))

    elseif prefix == 'U'
        push!(path, Point(current_ending.x, current_ending.y - distance))
    else
        error("Somehow, prefix broken")
    end
    return path
end

function orientation(inst::Tuple{Point, Point})
    if inst[1].x == inst[2].x  # x stays the same
        return 1
    else
        return 2
    end
end

function colinear(inst1::Tuple{Point, Point}, inst2::Tuple{Point, Point})
  #println(inst1, inst2)
  (x_dist, y_dist) = (Inf, Inf)
    if inst1[1].x >= min(inst2[1].x, inst2[2].x) && inst1[1].x <= max(inst2[1].x, inst2[2].x)
      x_dist = inst1[1].x
    elseif inst2[1].x >= min(inst1[1].x, inst1[2].x) && inst2[1].x <= max(inst1[1].x, inst1[2].x)
      x_dist = inst2[1].x
    end
    if inst1[1].y >= min(inst2[1].y, inst2[2].y) && inst1[1].y <= max(inst2[1].y, inst2[2].y)
      y_dist = inst1[1].y
    elseif inst2[1].y >= min(inst1[1].y, inst1[2].y) && inst2[1].y <= max(inst1[1].y, inst1[2].y)
      y_dist = inst2[1].y
    end
    #println("colinear", (x_dist, y_dist))
    return (x_dist, y_dist)
end

function cross_intersect(inst1::Tuple{Point, Point}, inst2::Tuple{Point, Point})
  # assume that inst1 has static x value and inst2 has static y value
  (x_dist, y_dist) = (Inf, Inf)
  #println(inst1, inst2)
  if inst1[1].x <= max(inst2[1].x, inst2[2].x) && inst1[1].x >= min(inst2[1].x, inst2[2].x) &&
    inst2[1].y <= max(inst1[1].y, inst1[2].y) && inst2[1].y >= min(inst1[1].y, inst1[2].y)
    x_dist = inst1[1].x
    y_dist = inst2[1].y
  end

  #println("cross", (x_dist, y_dist))
  return (x_dist, y_dist)
end

function test_intersect(inst1::Tuple{Point, Point}, inst2::Tuple{Point, Point}) :: Tuple{Bool, Point}
  x_dist = Inf
  y_dist = Inf
    if orientation(inst1) == orientation(inst2)
      (x_dist, y_dist) = colinear(inst1, inst2)
    else #different orientation
      #println(orientation(inst1))
      if orientation(inst1) == 1
        (x_dist, y_dist) = cross_intersect(inst1, inst2)
      else
        (x_dist, y_dist) = cross_intersect(inst2, inst1)
      end
    end

    if x_dist != Inf && y_dist != Inf
      return (true, Point(x_dist, y_dist))
    end
    return (false, Point(0,0))
end

function algorithm(input::Array{Array{String,1},1})
    inp1 = foldl(parse_prefix, input[1]; init=[Point(0,0)] )
    inp2 = foldl(parse_prefix, input[2]; init=[Point(0,0)] )
    best_guess = Inf
    for i = 2:length(inp1)-1
        for j = 2:length(inp2)-1
            entry1 = (inp1[i], inp1[i+1])
            entry2 = (inp2[j], inp2[j+1])
            check = test_intersect(entry1, entry2)
            #println(check)
            if check[1] == true
              best_guess = min(best_guess, manhattan_dist(check[2]))
            end
        end
    end
    best_guess
end

function manhattan_dist(p1::Point, p2::Point = Point(0,0)) :: Int64
  return abs(p2.x - p1.x) + abs(p2.y - p1.y)
end

function algorithm2(input::Array{Array{String,1},1}) :: Int64
    inp1 = foldl(parse_prefix, input[1]; init=[Point(0,0)] )
    inp2 = foldl(parse_prefix, input[2]; init=[Point(0,0)] )
    best_guess = Inf
    manhattan1 = [manhattan_dist(inp1[2], inp1[1])]
    for i = 2:length(inp1)-1
      manhattan2 = [manhattan_dist(inp2[2], inp2[1])]
        for j = 2:length(inp2)-1
            entry1 = (inp1[i], inp1[i+1])
            entry2 = (inp2[j], inp2[j+1])
            check = test_intersect(entry1, entry2)
            #println(check)
            if check[1] == true
              println(manhattan1)
              println(manhattan2)
              intersect_diff = manhattan_dist(inp1[i], check[2]) + manhattan_dist(inp2[j], check[2])
              #best_guess = min(best_guess, manhattan_dist(check[2]))
              best_guess = min(best_guess, sum(manhattan1) + sum(manhattan2) + intersect_diff)
            end
          push!(manhattan2, manhattan_dist(inp2[j+1], inp2[j]))
        end
      push!(manhattan1, manhattan_dist(inp1[i+1], inp1[i]))
    end
    best_guess
end

using Test
example_input1 = [["R8","U5","L5","D3"], ["U7","R6","D4","L4"]]
example_input2 = [["R75","D30","R83","U83","L12","D49","R71","U7","L72"],
["U62","R66","U55","R34","D71","R55","D58","R83"]]
@test parse_prefix([Point(0,0)], "R75") == [Point(0,0), Point(75,0)]
@test parse_prefix([Point(75, 0)], "D30") == [Point(75,0), Point(75, 30)]
@test algorithm(example_input1) == 6
@test algorithm(example_input2) == 159
@test algorithm2(example_input1) == 30
println(input_data)
println(algorithm(input_data))
println(algorithm2(input_data))

max_stars_repo_name: buotex/aoc

max_issues_repo_name: buotex/aoc

max_forks_repo_name: buotex/aoc


text:
"""
    IrfEstimate

`struct` to hold impulse response functions for a [`VAREstimate`](@ref)

# Fields
- point: Point estimates of impulse responses
- upper: Upper bounds of IRF confidence intervals
- lower: Lower bounds of IRF confidence intervals

"""
struct IrfEstimate
    point::Array
    upper::Array
    lower::Array
end

"""

"""
function IRF(
    v::VAREstimate,
    h::Int;
    # boot::Bool = false,
    # reps::Int = 500,
    ci::Float64 = 0.95,
    orthogonal = true,
    cumulative = false,
)
    Y = v.Y
    p = v.lags
    constant = v.constant
    trend = v.trend
    B = v.B
    U = v.U
    Σ = v.ΣU

    K = size(B, 2)
    ϕ = simple_irf(B, K, p, h)
    if orthogonal
        Θ = orthogonal_irf(ϕ, Σ, K, h)
        point = permutedims(Θ, [1, 3, 2])
    else
        point = permutedims(ϕ, [1, 3, 2])
    end

    # if boot
    #     upper, lower =
    #         bootstrap_irf_ci(Y, p, constant, trend, B, U, h; reps = reps, ci = ci)
    # else
    #     upper = zeros(K, h + 1, K)
    #     lower = zeros(K, h + 1, K)
    # end

    upper = zeros(K, h + 1, K)
    lower = zeros(K, h + 1, K)

    # Cumulative impulse responses
    if cumulative
        # upper = cumsum(upper, dims = 2)
        # lower = cumsum(lower, dims = 2)
        point = cumsum(point, dims = 2)
    end

    IrfEstimate(point, upper, lower)
end

function simple_irf(B, K, p, h)
    nparam = size(B, 1)
    A = B[(nparam-K*p)+1:end, :]
    A = reshape(A', (K, K, p))
    phi = zeros(K, K, h + 1)
    phi[:, :, 1] = Matrix(1.0I, K, K)
    @views for i = 2:h+1
        for j = 1:min(i - 1, p)
            @inbounds phi[:, :, i] += phi[:, :, i-j] * A[:, :, j]
        end
    end
    return phi
end

function orthogonal_irf(phi, sigma, K, h)
    # Orthogonalize using lower Cholesky factorization
    P = cholesky(sigma).L
    theta = zeros(K, K, h + 1)
    theta[:, :, 1] = P
    @views for i = 1:(h+1)
        theta[:, :, i] = phi[:, :, i] * P
    end
    return theta
end

# function bootstrap_irf_ci(Y, p, c, t, B, U, h; reps = 500, ci = 0.95)
#     nobs, K = size(Y)
#     uobs = nobs - p
#     ps = Y[1:p, :]
#     lb = (1 - ci) / 2
#     ub = 1 - lb
#     sim_sirf = []
#     sim_oirf = []
#     sim_irf_store = zeros(K * reps, h + 1, K)
#     upper = zeros(K, h + 1, K)
#     lower = zeros(K, h + 1, K)
#
#     for rep = 1:reps
#         Ysim = simulate_var(B, K, p, c, t, uobs; presample = ps, resid = U)
#         simVAR = VAR(Ysim, p; constant = c, trend = t)
#         simB = simVAR.B
#         simSigma = simVAR.ΣU
#         sim_sirf = simple_irf(simB, K, p, h + 1)
#         sim_oirf = orthogonalize_irf(sim_sirf, simSigma, K, h)
#         #=============
#         Store results:
#         -- "Layer" k (3rd dim) contains all responses to a shock to varible k
#         -- Within each layer, all replications for the response of each variable are grouped
#         together; e.g. rows [1:reps] contain all replications for the response of variable 1,
#         rows [reps+1:2*reps] contain all replications for the response of variable 2, etc.
#         ==============#
#         for sh = 1:K
#             for res = 1:K
#                 sim_irf_store[reps*(res-1)+rep, :, sh] = sim_oirf[res, sh, :]
#             end
#         end # storage loop
#
#         # Get quantiles
#         for sh = 1:K
#             for res = 1:K
#                 for hor = 1:h+1
#                     lower[res, hor, sh] =
#                         quantile(sim_irf_store[reps*(res-1)+1:(reps*res), hor, sh], lb)
#                     upper[res, hor, sh] =
#                         quantile(sim_irf_store[reps*(res-1)+1:(reps*res), hor, sh], ub)
#                 end
#             end
#         end # quantile loop
#     end # rep loop
#
#     return upper, lower
# end
#
# function simulate_var(
#     B,
#     K,
#     p,
#     constant,
#     trend,
#     uobs;
#     presample = nothing,
#     sigma = nothing,
#     resid = nothing,
# )
#     # Initialize simulated data
#     simY = zeros(uobs + p, K)
#     # Set presample values
#     if !isa(presample, Nothing)
#         simY[1:p, :] = presample
#     end
#     # Set covariance matrix for disturbances
#     if isa(sigma, Nothing)
#         Σ = Matrix(1.0I, K, K)
#     else
#         Σ = sigma
#     end
#     # Draw sample of residuals (with replacement) if provided
#     if !isa(resid, Nothing)
#         boot = true
#         randidx = sample(collect(1:uobs), uobs)
#         bootU = resid[randidx, :]
#     else
#         boot = false
#     end
#
#     for t = (p+1):(uobs+p)
#         # lagged endogenous variables for RHS
#         Zt = transpose(simY[t-1:-1:t-p, :])[:]
#         if trend
#             Zt = [t; Zt]
#         end
#         if constant
#             Zt = [1; Zt]
#         end
#
#         if boot
#             simY[t, :] = B' * Zt + bootU[t-p, :]
#         else
#             simY[t, :] = B' * Zt + rand(MvNormal(zeros(K), Σ))
#         end
#     end # rep loop
#
#     return simY
# end


text:
using MAT
using NPZ
using Dates
using GLPK
using JuMP
using Plots
using Plots.PlotMeasures
include("../utils.jl")

function get_wave_t_in_range(param, min, max)
    indices_to_keep = findall(>=(min), param)
    param = param[indices_to_keep]

    indices_to_keep = findall(<(max), param)
    param = param[indices_to_keep]
    param, indices_to_keep
end

dir_path = get_dataset_path(;dataset_name="grand_popo")
kalman_bathy_vars = matread(string(dir_path, "Kalman_bathy.mat"))
bathy_xb = kalman_bathy_vars["xb"][1, :]
bathy_tn = kalman_bathy_vars["tn"][1, :]
bathy_hbK = kalman_bathy_vars["hbK"]

min_target_matdate = minimum(filter(!isnan, kalman_bathy_vars["tn"]))
max_target_matdate = maximum(filter(!isnan, kalman_bathy_vars["tn"]))

wave_vars = matread(string(dir_path, "Wave_GPP_ERA.mat"))
wave_t = wave_vars["t"]
wave_Tm = wave_vars["Tm"]
wave_hs = wave_vars["hs"]
wave_dir = wave_vars["dir"]

wave_t, indices = get_wave_t_in_range(wave_t, min_target_matdate, max_target_matdate)
wave_Tm = wave_Tm[indices]
wave_hs = wave_hs[indices]
wave_dir = wave_dir[indices]

outdir = "/Users/mn/JuliaProjects/cgp-morphcast/output/grand_popo/data/"

npzwrite(string(outdir, "bathy_xb.npy"), bathy_xb)
npzwrite(string(outdir, "bathy_tn.npy"), bathy_tn)
npzwrite(string(outdir, "bathy_hbK.npy"), bathy_hbK)

npzwrite(string(outdir, "wave_t.npy"), wave_t)
npzwrite(string(outdir, "wave_Tm.npy"), wave_Tm)
npzwrite(string(outdir, "wave_hs.npy"), wave_hs)
npzwrite(string(outdir, "wave_dir.npy"), wave_dir)

max_stars_repo_name: mahmoud-al-najar/cgp-morphcast

max_issues_repo_name: mahmoud-al-najar/cgp-morphcast

max_forks_repo_name: mahmoud-al-najar/cgp-morphcast


text:
mutable struct FTestResult{N}
    ssr::NTuple{N, Float64}
    dof::NTuple{N, Int}
    dof_resid::NTuple{N, Int}
    r2::NTuple{N, Float64}
    fstat::Tuple{Vararg{Float64}}
    pval::Tuple{Vararg{PValue}}
end

"""A helper function to determine if mod1 is nested in mod2"""
function issubmodel(mod1::LinPredModel, mod2::LinPredModel; atol=0::Real)
    mod1.rr.y != mod2.rr.y && return false # Response variables must be equal

    # Test that models are nested
    pred1 = mod1.pp.X
    npreds1 = size(pred1, 2)
    pred2 = mod2.pp.X
    npreds2 = size(pred2, 2)
    # If model 1 has more predictors, it can't possibly be a submodel
    npreds1 > npreds2 && return false
    # Test min norm pred2*B - pred1 ≈ 0
    rtol = Base.rtoldefault(typeof(pred1[1,1]))
    nresp = size(pred2, 1)
    return norm(view(qr(pred2).Q'pred1, npreds2 + 1:nresp, :)) <= max(atol, rtol*norm(pred1))
end

_diffn(t::NTuple{N, T}) where {N, T} = ntuple(i->t[i]-t[i+1], N-1)

_diff(t::NTuple{N, T}) where {N, T} = ntuple(i->t[i+1]-t[i], N-1)

dividetuples(t1::NTuple{N}, t2::NTuple{N}) where {N} = ntuple(i->t1[i]/t2[i], N)

"""
    ftest(mod::LinearModel...; atol=0::Real)

For each sequential pair of linear predictors in `mod`, perform an F-test to determine if
the first one fits significantly better than the next.

A table is returned containing residual degrees of freedom (DOF), degrees of freedom,
difference in DOF from the preceding model, sum of squared residuals (SSR), difference in
SSR from the preceding model, R², difference in R² from the preceding model, and F-statistic
and p-value for the comparison between the two models.

!!! note
    This function can be used to perform an ANOVA by testing the relative fit of two models
    to the data

Optional keyword argument `atol` controls the numerical tolerance when testing whether
the models are nested.

# Examples

Suppose we want to compare the effects of two or more treatments on some result. Because
this is an ANOVA, our null hypothesis is that `Result ~ 1` fits the data as well as
`Result ~ 1 + Treatment`.

```jldoctest
julia> dat = DataFrame(Treatment=[1, 1, 1, 2, 2, 2, 1, 1, 1, 2, 2, 2.],
                       Result=[1.1, 1.2, 1, 2.2, 1.9, 2, .9, 1, 1, 2.2, 2, 2],
                       Other=[1, 1, 2, 1, 2, 1, 3, 1, 1, 2, 2, 1]);

julia> model = lm(@formula(Result ~ 1 + Treatment), dat);

julia> nullmodel = lm(@formula(Result ~ 1), dat);

julia> bigmodel = lm(@formula(Result ~ 1 + Treatment + Other), dat);

julia> ft = ftest(model.model, nullmodel.model)
        Res. DOF DOF ΔDOF    SSR    ΔSSR     R²    ΔR²       F* p(>F)
Model 1       10   3      0.1283         0.9603
Model 2       11   2   -1 3.2292 -3.1008 0.0000 0.9603 241.6234 <1e-7

julia> ftest(bigmodel.model, model.model, nullmodel.model)
        Res. DOF DOF ΔDOF    SSR    ΔSSR     R²    ΔR²       F*  p(>F)
Model 1        9   4      0.1038         0.9678
Model 2       10   3   -1 0.1283 -0.0245 0.9603 0.0076   2.1236 0.1790
Model 3       11   2   -1 3.2292 -3.1008 0.0000 0.9603 241.6234  <1e-7
```
"""
function ftest(mods::LinearModel...; atol=0::Real)
    nmodels = length(mods)
    for i in 2:nmodels
        issubmodel(mods[i], mods[i-1], atol=atol) ||
        throw(ArgumentError("F test $i is only valid if model $i is nested in model $(i-1)"))
    end

    SSR = deviance.(mods)

    nparams = dof.(mods)

    df2 = _diffn(nparams)
    df1 = Int.(dof_residual.(mods))

    MSR1 = dividetuples(_diff(SSR), df2)
    MSR2 = dividetuples(SSR, df1)[1:nmodels-1]

    fstat = dividetuples(MSR1, MSR2)
    pval = PValue.(ccdf.(FDist.(df2, df1[1:nmodels-1]), fstat))
    return FTestResult(SSR, nparams, df1, r2.(mods), fstat, pval)
end

function show(io::IO, ftr::FTestResult{N}) where N
    Δdof = _diffn(ftr.dof_resid)
    Δssr = _diffn(ftr.ssr)
    ΔR² = _diffn(ftr.r2)

    nc = 10
    nr = N
    outrows = Matrix{String}(undef, nr+1, nc)

    outrows[1, :] = ["", "Res. DOF", "DOF", "ΔDOF", "SSR", "ΔSSR",
                     "R²", "ΔR²", "F*", "p(>F)"]

    outrows[2, :] = ["Model 1", @sprintf("%.0d", ftr.dof_resid[1]),
                     @sprintf("%.0d", ftr.dof[1]), " ",
                     @sprintf("%.4f", ftr.ssr[1]), " ",
                     @sprintf("%.4f", ftr.r2[1]), " ", " ", " "]

    for i in 2:nr
        outrows[i+1, :] = ["Model $i", @sprintf("%.0d", ftr.dof_resid[i]),
                           @sprintf("%.0d", ftr.dof[i]), @sprintf("%.0d", Δdof[i-1]),
                           @sprintf("%.4f", ftr.ssr[i]), @sprintf("%.4f", Δssr[i-1]),
                           @sprintf("%.4f", ftr.r2[i]), @sprintf("%.4f", ΔR²[i-1]),
                           @sprintf("%.4f", ftr.fstat[i-1]), string(ftr.pval[i-1]) ]
    end
    colwidths = length.(outrows)
    max_colwidths = [maximum(view(colwidths, :, i)) for i in 1:nc]

    for r in 1:nr+1
        for c in 1:nc
            cur_cell = outrows[r, c]
            cur_cell_len = length(cur_cell)

            padding = " "^(max_colwidths[c]-cur_cell_len)
            if c > 1
                padding = " "*padding
            end

            print(io, padding)
            print(io, cur_cell)
        end
        print(io, "\n")
    end
end


text:
using BenchmarkTools
using Printf
using Transducers

const BoardingPass = NamedTuple{(:row, :column, :seatId), Tuple{Int64, Int64, Int64}}

function toInt(s::AbstractString, zeroChar::Char, oneChar::Char)::Int64
    value = 0
    for i = 1:length(s)
        value *= 2
        if s[i] == oneChar
            value += 1
        elseif s[i] != zeroChar
            throw(DomainError(s, "valid characters are $zeroChar and $oneChar"))
        end
    end
    value
end

function toBoardingPass(s::AbstractString)::BoardingPass
    row = toInt(s[1:7], 'F', 'B')
    column = toInt(s[8:10], 'L', 'R')
    (row = row, column = column, seatId = row*8+column)
end

function toBoardingPass(seatId::Int64)::BoardingPass
    (row = seatId ÷ 8, column = seatId % 8, seatId = seatId)
end

#### basic tests

@assert toBoardingPass("FBFBBFFRLR") == (row = 44, column = 5, seatId = 357)
@assert toBoardingPass("BFFFBBFRRR") == (row = 70, column = 7, seatId = 567)
@assert toBoardingPass("FFFBBBFRRR") == (row = 14, column = 7, seatId = 119)
@assert toBoardingPass("BBFFBBFRLL") == (row = 102, column = 4, seatId = 820)

@assert toBoardingPass(357) == (row = 44, column = 5, seatId = 357)
@assert toBoardingPass(567) == (row = 70, column = 7, seatId = 567)
@assert toBoardingPass(119) == (row = 14, column = 7, seatId = 119)
@assert toBoardingPass(820) == (row = 102, column = 4, seatId = 820)

#### read data

boardingPasses = readlines("src/advent/day5/input.txt") |> Map(line -> toBoardingPass(line)) |> collect

#### part 1

function highestSeatId(boardingPasses)
    boardingPasses |> Map(pass -> pass.seatId) |> maximum
end

@printf("highestSeatId = %d\n", highestSeatId(boardingPasses))
#@btime highestSeatId(boardingPasses);

#### part 2

function findMySeatId(boardingPasses)
    map = Dict{Int64, BoardingPass}()
    for pass in boardingPasses
        map[pass.seatId] = pass
    end
    for pass in boardingPasses
        if haskey(map, pass.seatId - 2) && !haskey(map, pass.seatId - 1)
            return toBoardingPass(pass.seatId - 1)
        elseif haskey(map, pass.seatId + 2) && !haskey(map, pass.seatId + 1)
            return toBoardingPass(pass.seatId + 1)
        end
    end
end

@printf("mySeatId = %d\n", findMySeatId(boardingPasses).seatId)
#@btime findMySeatId(boardingPasses);

#### part 2: pure MapReduce

function merge(a::UnitRange{Int64}, b::UnitRange{Int64})::Tuple{UnitRange{Int64}, Vararg{UnitRange{Int64}}}
    a.stop + 1 == b.start && return (a.start:b.stop, )
    b.stop + 1 == a.start && return (b.start:a.stop, )
    (a, b)
end

@assert merge(1:1, 2:2) == (1:2, )
@assert merge(2:2, 1:1) == (1:2, )
@assert merge(2:2, 4:4) == (2:2, 4:4)

function merge!(ranges, range1)
    modified = false
    for i ∈ 1:length(ranges)
        results = merge(ranges[i], range1)
        if length(results) == 1
            ranges[i] = results[1]
            modified = true
            break
        end
    end
    if !modified
        push!(ranges, range1)
    end
    ranges
end

function merge!(ranges)
    modified = false
    while length(ranges) > 1
        changed = false
        i = 1
        while i < length(ranges)
            range = ranges[i]
            for j = i+1:length(ranges)
                results = merge(range, ranges[j])
                if length(results) == 1
                    ranges[i] = results[1]
                    deleteat!(ranges, j)
                    changed = true
                    modified = true
                    break
                end
            end
            i += 1
        end
        !changed && break
    end
    modified
end

function mergeRangeList(a, b)
    #before = (length(a), length(b))
    for range ∈ b
        merge!(a, range)
    end
    while merge!(a)
        # empty body
    end
    #after = length(a)
    #println("---- $before -> $after")
    a
end

function findMySeatId2(boardingPasses)
    ranges = foldxt(mergeRangeList, boardingPasses |> Map(pass -> [pass.seatId:pass.seatId])) |> collect
    sort!(ranges)
    1:length(ranges)-1 |> Map(i -> (ranges[i], ranges[i + 1])) |> Filter(ranges -> ranges[2].start - ranges[1].stop == 2) |> Map(ranges -> toBoardingPass(ranges[1].stop + 1)) |> first

    #for i = 1:length(ranges)-1
    #    range1 = ranges[i]
    #    range2 = ranges[i + 1]
    #    range2.start - range1.stop == 2 && return (seatId = range1.stop + 1, found = true)
    #end
    #(seatId = -1, found = false)
end

@printf("mySeatId2 = %d\n", findMySeatId2(boardingPasses).seatId)
#@btime findMySeatId2(boardingPasses);


text:
"""
Used to provide a definition for a Gaussian mutation operator.
"""
type GaussianMutationDefinition <: MutationDefinition
  stage::AbstractString
  rate::Float
  mean::Float
  std::Float

  GaussianMutationDefinition() = new("", 0.01, 0.0, 1.0)
  GaussianMutationDefinition(stage::AbstractString) = new(stage, 0.01, 0.0, 1.0)
  GaussianMutationDefinition(r::Float) = new("", r, 0.0, 1.0)
  GaussianMutatianDefinition(r::Float, u::Float, s::Float) =
    new("", r, u, s)
  GaussianMutationDefinition(stage::AbstractString, r::Float, u::Float, s::Float) =
    new(stage, r, u, s)
end

"""
Composes a Gaussian mutation operator from a provided definition.
"""
function compose!(d::GaussianMutationDefinition, sp::Species)
  d.stage = d.stage == "" ? genotype(sp).label : d.stage
  GaussianMutation{Float}(d.mean, d.std, d.rate, r)
end

type GaussianMutation{T <: Number} <: Mutation
  stage::AbstractString
  mean::T
  std::T
  min_value::T
  max_value::T
  rate::Float
  representation::Representation

  GaussianMutation(stage::AbstractString, mean::T, std::T, rate::Float, rep::Representation) =
    new(stage, mean, std, representation.minimum_value(rep), representation.maximum_value(rep), rate, rep)
end

"""
TODO: DESCRIPTION OF GAUSSIAN MUTATION.
"""
gaussian() = GaussianMutationDefinition()
gaussian(rate::Float) = GaussianMutationDefinition(rate)
gaussian(stage::AbstractString) = GaussianMutationDefinition(stage)
gaussian(rate::Float, mean::Float, std::Float) =
  GaussianMutationDefinition(rate, mean, std)
gaussian(stage::AbstractString, rate::Float, mean::Float, std::Float) =
  GaussianMutationDefinition(stage, rate, mean, std)
function gaussian(f::Function)
  def = GaussianMutationDefinition()
  f(def)
  def
end

"""
Produces some Gaussian noise from a given distribution, specified by its mean
and standard deviation.
"""
gaussian_noise{T <: Number}(mean::T, std::T) =
  mean + std * sqrt(-2 * log(1 - rand())) * sin(2 * pi * rand())

function mutate!{T <: Number}(o::GaussianMutation, input::Vector{T})
  # If we knew the length of our genome in advance, how much faster
  # would this operation be?
  for i in 1:length(input)
    if rand() <= o.rate
      input[i] += gaussian_noise(o.mean, o.std)
      input[i] = min(max(input[i], o.min_value), o.max_value)
    end
  end
  return input
end


text:
module HERAEUS
using ..AGFFileReader: Glass, GlassID, AGF
export HERASIL, HOMOSIL, SUPRASIL, INFRASIL, HOQ

const HERASIL = Glass(GlassID(AGF, 1254), 2, 0.47652307, 0.00284888095, 0.627786368, 0.0118369052, 0.872274404, 95.6856012, 0.0, 0.0, NaN, NaN, 0.20254, 3.5, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 20.0, 0.0, -1.0, -1.0, 0.51, -1.0, 0.0, -1.0, [], 1.458563, -1.0, -1.0, 0, 67.664757, 0, 2.203, 0.0)
const HOMOSIL = Glass(GlassID(AGF, 1255), 2, 0.47652307, 0.00284888095, 0.627786368, 0.0118369052, 0.872274404, 95.6856012, 0.0, 0.0, NaN, NaN, 0.20254, 3.5, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 20.0, 0.0, -1.0, -1.0, 0.51, -1.0, 0.0, -1.0, [], 1.458563, -1.0, -1.0, 0, 67.664757, 0, 2.203, 0.0)
const SUPRASIL = Glass(GlassID(AGF, 1256), 2, 0.473115591, 0.012995717, 0.631038719, 0.0041280922, 0.906404498, 98.7685322, 0.0, 0.0, NaN, NaN, 0.19, 3.5, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 20.0, 0.0, -1.0, -1.0, 0.51, -1.0, 0.0, -1.0, [], 1.458464, -1.0, -1.0, 0, 67.878205, 0, 2.201, 0.0)
const INFRASIL = Glass(GlassID(AGF, 1257), 2, 0.47652307, 0.00284888095, 0.627786368, 0.0118369052, 0.872274404, 95.6856012, 0.0, 0.0, NaN, NaN, 0.20254, 3.5, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 20.0, 0.0, -1.0, -1.0, 0.51, -1.0, 0.0, -1.0, [], 1.458563, -1.0, -1.0, 0, 67.664757, 0, 2.203, 0.0)
const HOQ = Glass(GlassID(AGF, 1258), 2, 0.47652307, 0.00284888095, 0.627786368, 0.0118369052, 0.872274404, 95.6856012, 0.0, 0.0, NaN, NaN, 0.20254, 3.5, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 20.0, 0.0, -1.0, -1.0, 0.51, -1.0, 0.0, -1.0, [], 1.458563, -1.0, -1.0, 0, 67.664757, 0, 2.203, 0.0)
end # module

max_stars_repo_name: rambunctiousapple/GlassCat

max_issues_repo_name: rambunctiousapple/GlassCat

max_forks_repo_name: rambunctiousapple/GlassCat


text:
function equals(obj::SimplexTableau, arg0::Object)
    return jcall(obj, "equals", jboolean, (Object,), arg0)
end

function hash_code(obj::SimplexTableau)
    return jcall(obj, "hashCode", jint, ())
end

function normalize_constraints(obj::SimplexTableau, arg0::Collection)
    return jcall(obj, "normalizeConstraints", List, (Collection,), arg0)
end


text:
function bisect(X::Interval)
    m = mid(X)
    return ( Interval(X.lo, m), Interval(m, X.hi) )
end

function bisect(X::IntervalBox)
    i = findmax([diam(x) for x in X])[2]  # find largest

    x1, x2 = bisect(X[i])

    # the following is ugly -- replace by iterators once https://github.com/JuliaLang/julia/pull/15516
    # has landed?
    # X1 = IntervalBox(X[1:i-1]..., x1, X[i+1:end]...)  # insert x1 in i'th place
    # X2 = IntervalBox(X[1:i-1]..., x2, X[i+1:end]...)
    X1 = setindex(X, x1, i)
    X2 = setindex(X, x2, i)

    return (X1, X2)
end

doc"""
`pave` takes the given working list of boxes and splits them into inner and boundary
lists with the given separator
"""
function pave{N,T}(S::Separator, working::Vector{IntervalBox{N,T}}, ϵ)

    inner_list = SubPaving{N,T}()
    boundary_list = SubPaving{N,T}()

    while !isempty(working)

        X = pop!(working)

        inner, outer = S(X)   # here inner and outer are reversed compared to Jaulin
        # S(X) returns the pair (contractor with respect to the inside of the constraing, contractor with respect to outside)

        inner2 = IntervalBox(inner)
        outer2 = IntervalBox(outer)

        inside_list = setdiff(X, outer2)

        if length(inside_list) > 0
            append!(inner_list, inside_list)
        end

        boundary = inner2 ∩ outer2

        if isempty(boundary)
            continue
        end

        if diam(boundary) < ϵ
            push!(boundary_list, boundary)

        else
            push!(working, bisect(boundary)...)
        end

    end

    return inner_list, boundary_list

end

doc"""
    pave(S::Separator, domain::IntervalBox, eps)`

Find the subset of `domain` defined by the constraints specified by the separator `S`.
Returns (sub)pavings `inner` and `boundary`, i.e. lists of `IntervalBox`.
"""
function pave{N,T}(S::Separator, X::IntervalBox{N,T}, ϵ = 1e-2)

    inner_list, boundary_list = pave(S, [X], ϵ)

    return Paving{N,T}(S, inner_list, boundary_list, ϵ)

end

doc"""Refine a paving to tolerance ϵ"""
function refine!(P::Paving, ϵ = 1e-2)
    if P.ϵ <= ϵ  # already refined
        return
    end

    new_inner, new_boundary = pave(P.separator, P.boundary, ϵ)

    append!(P.inner, new_inner)
    P.boundary = new_boundary
    P.ϵ = ϵ
end


text:
function sol_min(result, df_result)

    for bus_key in [i for i in keys(result["solution"]["bus"])]

        va_key = "va" * "_" * bus_key
        va_value = []
        push!(va_value, string(result["solution"]["bus"][bus_key]["va"]))

        vm_key = "vm" * "_" * bus_key
        vm_value = []
        push!(vm_value, string(result["solution"]["bus"][bus_key]["vm"]))

        df_result[!, va_key] = va_value
        df_result[!, vm_key] = vm_value

    end

    for gen_key in [i for i in keys(result["solution"]["gen"])]

        qg_key = "qg" * "_" * gen_key
        qg_value = []
        push!(qg_value, string(result["solution"]["gen"][gen_key]["qg"]))

        pg_key = "pg" * "_" * gen_key
        pg_value = []
        push!(pg_value, string(result["solution"]["gen"][gen_key]["pg"]))

        df_result[!, qg_key] = qg_value
        df_result[!, pg_key] = pg_value

    end

    for branch_key in [i for i in keys(result["solution"]["branch"])]
        qf_key = "qf" * "_" * branch_key
        qf_value = []
        push!(qf_value, string(result["solution"]["branch"][branch_key]["qf"]))

        qt_key = "qt" * "_" * branch_key
        qt_value = []
        push!(qt_value, string(result["solution"]["branch"][branch_key]["qt"]))

        pf_key = "pf" * "_" * branch_key
        pf_value = []
        push!(pf_value, string(result["solution"]["branch"][branch_key]["pf"]))

        pt_key = "pt" * "_" * branch_key
        pt_value = []
        push!(pt_value, string(result["solution"]["branch"][branch_key]["pt"]))

        df_result[!, qf_key] = qf_value
        df_result[!, qt_key] = qt_value
        df_result[!, pf_key] = pf_value
        df_result[!, pt_key] = pt_value

    end

end

max_stars_repo_name: zhouliang-yu/Grid-ML

max_issues_repo_name: zhouliang-yu/Grid-ML

max_forks_repo_name: zhouliang-yu/Grid-ML


text:
function irepeat(v::AbstractVector, n::Int)
    nV = length(v)
    res = similar(v, nV * n)
    @inbounds for j = 1:nV
        vj = v[j]
        base = (j - 1) * n
        @inbounds @simd for i = 1:n
            res[base+i] = vj
        end
    end
    res
end

function orepeat(v::AbstractVector, n::Int)
    nV = length(v)
    res = similar(v, nV * n)
    @inbounds for i = 1:n
        base = (i - 1) * nV
        @inbounds @simd for j = 1:nV
            res[base+j] = v[j]
        end
    end
    res
end

# TODO: since 0.7 transpose is different, we don't take transpose serious here.
####### kronecker product ###########
# TODO: if IMatrix{1}, do nothing
kron(A::IMatrix{Na,Ta}, B::IMatrix{Nb,Tb}) where {Na,Nb,Ta<:Number,Tb<:Number} =
    IMatrix{Na * Nb,promote_type(Ta, Tb)}()
kron(A::IMatrix{Na,Ta}, B::IMatrix{1,Tb}) where {Na,Ta<:Number,Tb<:Number} =
    IMatrix{Na,promote_type(Ta, Tb)}()
kron(A::IMatrix{1,Ta}, B::IMatrix{Nb,Tb}) where {Nb,Ta<:Number,Tb<:Number} =
    IMatrix{Nb,promote_type(Ta, Tb)}()
kron(A::IMatrix{Na,<:Number}, B::Diagonal{<:Number}) where {Na} =
    Diagonal(orepeat(B.diag, Na))
kron(B::Diagonal{<:Number}, A::IMatrix{Na}) where {Na} = Diagonal(irepeat(B.diag, Na))
for MT in [:AbstractMatrix, :PermMatrix, :SparseMatrixCSC, :Diagonal]
    @eval kron(A::IMatrix{1,<:Number}, B::$MT{<:Number}) = B
    @eval kron(B::$MT{<:Number}, A::IMatrix{1,<:Number}) = B
end

####### diagonal kron ########
kron(A::StridedMatrix{<:Number}, B::Diagonal{<:Number}) = kron(A, PermMatrix(B))
kron(A::Diagonal{<:Number}, B::StridedMatrix{<:Number}) = kron(PermMatrix(A), B)
kron(A::Diagonal{<:Number}, B::SparseMatrixCSC{<:Number}) = kron(PermMatrix(A), B)
kron(A::SparseMatrixCSC{<:Number}, B::Diagonal{<:Number}) = kron(A, PermMatrix(B))

function kron(A::AbstractMatrix{Tv}, B::IMatrix{Nb}) where {Nb,Tv<:Number}
    mA, nA = size(A)
    nzval = Vector{Tv}(undef, Nb * mA * nA)
    rowval = Vector{Int}(undef, Nb * mA * nA)
    colptr = collect(1:mA:Nb*mA*nA+1)
    @inbounds for j = 1:nA
        source = view(A, :, j)
        startbase = (j - 1) * Nb * mA - mA
        for j2 = 1:Nb
            start = startbase + j2 * mA
            row = j2 - Nb
            @inbounds @simd for i = 1:mA
                nzval[start+i] = source[i]
                rowval[start+i] = row + Nb * i
            end
        end
    end
    SparseMatrixCSC(mA * Nb, nA * Nb, colptr, rowval, nzval)
end

function kron(A::IMatrix{Na}, B::AbstractMatrix{Tv}) where {Na,Tv<:Number}
    mB, nB = size(B)
    rowval = Vector{Int}(undef, nB * mB * Na)
    nzval = Vector{Tv}(undef, nB * mB * Na)
    @inbounds for j = 1:Na
        r0 = (j - 1) * mB
        for j2 = 1:nB
            start = ((j - 1) * nB + j2 - 1) * mB
            @inbounds @simd for i = 1:mB
                rowval[start+i] = r0 + i
                nzval[start+i] = B[i, j2]
            end
        end
    end
    colptr = collect(1:mB:nB*mB*Na+1)
    SparseMatrixCSC(mB * Na, Na * nB, colptr, rowval, nzval)
end

function kron(A::IMatrix{Na}, B::SparseMatrixCSC{T}) where {Na,T<:Number}
    mB, nB = size(B)
    nV = nnz(B)
    nzval = Vector{T}(undef, Na * nV)
    rowval = Vector{Int}(undef, Na * nV)
    colptr = Vector{Int}(undef, nB * Na + 1)
    nzval = Vector{T}(undef, Na * nV)
    colptr[1] = 1
    for i = 1:Na
        r0 = (i - 1) * mB
        start = nV * (i - 1)
        @inbounds @simd for k = 1:nV
            rowval[start+k] = B.rowval[k] + r0
            nzval[start+k] = B.nzval[k]
        end
        colbase = (i - 1) * nB
        @inbounds @simd for j = 2:nB+1
            colptr[colbase+j] = B.colptr[j] + start
        end
    end
    SparseMatrixCSC(mB * Na, nB * Na, colptr, rowval, nzval)
end

function kron(A::SparseMatrixCSC{T}, B::IMatrix{Nb}) where {T<:Number,Nb}
    mA, nA = size(A)
    nV = nnz(A)
    rowval = Vector{Int}(undef, Nb * nV)
    colptr = Vector{Int}(undef, nA * Nb + 1)
    nzval = Vector{T}(undef, Nb * nV)
    z = 1
    colptr[1] = 1
    @inbounds for i = 1:nA
        rstart = A.colptr[i]
        rend = A.colptr[i+1] - 1
        colbase = (i - 1) * Nb + 1
        @inbounds for k = 1:Nb
            irow_Nb = k - Nb
            @inbounds @simd for r = rstart:rend
                rowval[z] = A.rowval[r] * Nb + irow_Nb
                nzval[z] = A.nzval[r]
                z += 1
            end
            colptr[colbase+k] = z
        end
    end
    SparseMatrixCSC(mA * Nb, nA * Nb, colptr, rowval, nzval)
end

function kron(A::PermMatrix{T}, B::IMatrix) where {T<:Number}
    nA = size(A, 1)
    nB = size(B, 1)
    vals = Vector{T}(undef, nB * nA)
    perm = Vector{Int}(undef, nB * nA)
    @inbounds for i = 1:nA
        start = (i - 1) * nB
        permAi = (A.perm[i] - 1) * nB
        val = A.vals[i]
        @inbounds @simd for j = 1:nB
            perm[start+j] = permAi + j
            vals[start+j] = val
        end
    end
    PermMatrix(perm, vals)
end

function kron(A::IMatrix, B::PermMatrix{Tv,Ti}) where {Tv<:Number,Ti<:Integer}
    nA = size(A, 1)
    nB = size(B, 1)
    perm = Vector{Int}(undef, nB * nA)
    vals = Vector{Tv}(undef, nB * nA)
    @inbounds for i = 1:nA
        start = (i - 1) * nB
        @inbounds @simd for j = 1:nB
            perm[start+j] = start + B.perm[j]
            vals[start+j] = B.vals[j]
        end
    end
    PermMatrix(perm, vals)
end

function kron(A::StridedMatrix{Tv}, B::PermMatrix{Tb}) where {Tv<:Number,Tb<:Number}
    mA, nA = size(A)
    nB = size(B, 1)
    perm = fast_invperm(B.perm)
    nzval = Vector{promote_type(Tv, Tb)}(undef, mA * nA * nB)
    rowval = Vector{Int}(undef, mA * nA * nB)
    colptr = collect(1:mA:nA*nB*mA+1)
    z = 1
    @inbounds for j = 1:nA
        @inbounds for j2 = 1:nB
            p2 = perm[j2]
            val2 = B.vals[p2]
            ir = p2
            @inbounds @simd for i = 1:mA
                nzval[z] = A[i, j] * val2  # merge
                rowval[z] = ir
                z += 1
                ir += nB
            end
        end
    end
    SparseMatrixCSC(mA * nB, nA * nB, colptr, rowval, nzval)
end

function kron(A::PermMatrix{Ta}, B::StridedMatrix{Tb}) where {Tb<:Number,Ta<:Number}
    mB, nB = size(B)
    nA = size(A, 1)
    perm = fast_invperm(A.perm)
    nzval = Vector{promote_type(Ta, Tb)}(undef, mB * nA * nB)
    rowval = Vector{Int}(undef, mB * nA * nB)
    colptr = collect(1:mB:nA*nB*mB+1)
    z = 0
    @inbounds for j = 1:nA
        colbase = (j - 1) * nB
        p1 = perm[j]
        val2 = A.vals[p1]
        ir = (p1 - 1) * mB
        for j2 = 1:nB
            @inbounds @simd for i2 = 1:mB
                nzval[z+i2] = B[i2, j2] * val2
                rowval[z+i2] = ir + i2
            end
            z += mB
        end
    end
    SparseMatrixCSC(nA * mB, nA * nB, colptr, rowval, nzval)
end

function kron(A::PermMatrix{<:Number}, B::PermMatrix{<:Number})
    nA = size(A, 1)
    nB = size(B, 1)
    vals = kron(A.vals, B.vals)
    perm = Vector{Int}(undef, nB * nA)
    @inbounds for i = 1:nA
        start = (i - 1) * nB
        permAi = (A.perm[i] - 1) * nB
        @inbounds @simd for j = 1:nB
            perm[start+j] = permAi + B.perm[j]
        end
    end
    PermMatrix(perm, vals)
end

kron(A::PermMatrix{<:Number}, B::Diagonal{<:Number}) = kron(A, PermMatrix(B))
kron(A::Diagonal{<:Number}, B::PermMatrix{<:Number}) = kron(PermMatrix(A), B)

function kron(A::PermMatrix{Ta}, B::SparseMatrixCSC{Tb}) where {Ta<:Number,Tb<:Number}
    nA = size(A, 1)
    mB, nB = size(B)
    nV = nnz(B)
    perm = fast_invperm(A.perm)
    nzval = Vector{promote_type(Ta, Tb)}(undef, nA * nV)
    rowval = Vector{Int}(undef, nA * nV)
    colptr = Vector{Int}(undef, nA * nB + 1)
    colptr[1] = 1
    @inbounds @simd for i = 1:nA
        start_row = (i - 1) * nV
        start_ri = (perm[i] - 1) * mB
        v0 = A.vals[perm[i]]
        @inbounds @simd for j = 1:nV
            nzval[start_row+j] = B.nzval[j] * v0
            rowval[start_row+j] = B.rowval[j] + start_ri
        end
        start_col = (i - 1) * nB + 1
        start_ci = (i - 1) * nV
        @inbounds @simd for j = 1:nB
            colptr[start_col+j] = B.colptr[j+1] + start_ci
        end
    end
    SparseMatrixCSC(mB * nA, nB * nA, colptr, rowval, nzval)
end

function kron(A::SparseMatrixCSC{T}, B::PermMatrix{Tb}) where {T<:Number,Tb<:Number}
    nB = size(B, 1)
    mA, nA = size(A)
    nV = nnz(A)
    perm = fast_invperm(B.perm)
    rowval = Vector{Int}(undef, nB * nV)
    colptr = Vector{Int}(undef, nA * nB + 1)
    nzval = Vector{promote_type(T, Tb)}(undef, nB * nV)
    z = 1
    colptr[z] = 1
    @inbounds for i = 1:nA
        rstart = A.colptr[i]
        rend = A.colptr[i+1] - 1
        @inbounds for k = 1:nB
            irow = perm[k]
            bval = B.vals[irow]
            irow_nB = irow - nB
            @inbounds @simd for r = rstart:rend
                rowval[z] = A.rowval[r] * nB + irow_nB
                nzval[z] = A.nzval[r] * bval
                z += 1
            end
            colptr[(i-1)*nB+k+1] = z
        end
    end
    SparseMatrixCSC(mA * nB, nA * nB, colptr, rowval, nzval)
end


text:
# ---
# title: 1110. Delete Nodes And Return Forest
# id: problem1110
# author: Tian Jun
# date: 2020-10-31
# difficulty: Medium
# categories: Tree, Depth-first Search
# link: <https://leetcode.com/problems/delete-nodes-and-return-forest/description/>
# hidden: true
# ---
#
# Given the `root` of a binary tree, each node in the tree has a distinct value.
#
# After deleting all nodes with a value in `to_delete`, we are left with a
# forest (a disjoint union of trees).
#
# Return the roots of the trees in the remaining forest.  You may return the
# result in any order.
#
#
#
# **Example 1:**
#
# **![](https://assets.leetcode.com/uploads/2019/07/01/screen-
# shot-2019-07-01-at-53836-pm.png)**
#
#
#
#     Input: root = [1,2,3,4,5,6,7], to_delete = [3,5]
#     Output: [[1,2,null,4],[6],[7]]
#
#
#
#
# **Constraints:**
#
#   * The number of nodes in the given tree is at most `1000`.
#   * Each node has a distinct value between `1` and `1000`.
#   * `to_delete.length <= 1000`
#   * `to_delete` contains distinct values between `1` and `1000`.
#
#
## @lc code=start
using LeetCode

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


text:
mutable struct BufferedIRCurve{C<:AbstractIRCurve} <: AbstractIRCurve
    buffer::Dict{Date, Float64}
    curve::C
end

BufferedIRCurve(curve::AbstractIRCurve) = BufferedIRCurve(Dict{Date, Float64}(), curve)

for fun in (:curve_get_name, :curve_get_daycount, :curve_get_compounding, :curve_get_method, :curve_get_date, :curve_get_dtm, :curve_get_zero_rates, :curve_get_model_parameters, :curve_get_spline_fit_on_rates, :curve_get_spline_fit_on_discount_factors)
    @eval begin
        ($fun)(curve::BufferedIRCurve) = ($fun)(curve.curve)
    end
end

curve_get_dict_parameter(curve::BufferedIRCurve, sym::Symbol) = curve_get_dict_parameter(curve.curve, sym)

function zero_rate(curve::BufferedIRCurve, maturity::Date)
    if !haskey(curve.buffer, maturity)
        curve.buffer[maturity] = _zero_rate(curve_get_method(curve), curve, maturity)
    end
    return curve.buffer[maturity]
end

function zero_rate(curve::BufferedIRCurve, maturity_vec::Vector{Date})
    n = length(maturity_vec)
    result = Vector{Float64}(undef, n)
    for i in 1:n
        @inbounds result[i] = zero_rate(curve, maturity_vec[i])
    end
    return result
end


text:
# using Dates
using Base.Filesystem
needSaveFile=true
needJuliaSend=true

# println("000")
@static if Base.Sys.islinux()
    root=homedir()
    root = pwd()
    println( "root is :$root")
    Logfolder= joinpath( root,"../dataroot/","log/")

    folderA = joinpath( root,"../dataroot/","b/")
    folderB = joinpath( root,"../dataroot/","a/")
else
    root="d:\\"
    root = pwd()
    Logfolder= joinpath( root,"../dataroot\\","log\\")

    folderA = joinpath( root,"../dataroot\\","b\\")
    folderB = joinpath( root,"../dataroot\\","a\\")
end

# dataRoot=joinpath(root,"dataroot/")
println("Root data Folder is:",root )

IpStr = "192.168.99.193";
IpStr = "192.168.0.109"
IpStr = "127.0.0.1"
IpStr="10.14.38.51"
#IpStr = "172.20.160.131"
#IpStr = "192.168.1.5"
IpStr = "0.0.0.0"
#localhost
Port = UInt16(19998);
# MBytes
Unit= 1024*1024
# should be in Mib ,1G=124Mib
DiskTotalSpace = 80*1024

# 特征数据可占用总磁盘的百分比例，包含RAW3~4.其它的为 RAW1~2 占用
FeatureRatio=80

RAW12Total = DiskTotalSpace *(100-FeatureRatio ) /100
RAW34Total = DiskTotalSpace *FeatureRatio /100

RAW1MaxSpace= floor( UInt32,RAW12Total / 2 *2)
RAW2MaxSpace= floor( UInt32,RAW12Total / 2*0 )
RAW3MaxSpace= floor( UInt32,RAW34Total / 2 *2 )
RAW4MaxSpace= floor( UInt32,RAW34Total/  2 *0 )
# space = getDiskSpace( root )
###################################
for pa in (Logfolder, folderA, folderB )
    Filesystem.mkpath( pa )
end
struct cfgDataOnline
    folderA
    folderB
    IpStr
    Port
    Unit
    RAW1MaxSpace
    RAW2MaxSpace
    RAW3MaxSpace
    RAW4MaxSpace
    needSaveFile
    needJuliaSend
    Logfolder
    function cfgDataOnline(folderA,folderB ,IpStr,Port,Unit,RAW1MaxSpace,RAW2MaxSpace,RAW3MaxSpace,RAW4MaxSpace,needSaveFile,needJuliaSend,Logfolder)
        new(folderA,folderB ,IpStr,Port,Unit,RAW1MaxSpace,RAW2MaxSpace,RAW3MaxSpace,RAW4MaxSpace,needSaveFile,needJuliaSend,Logfolder)
    end
end

cfgObj = cfgDataOnline(folderA,folderB ,IpStr,Port,Unit,RAW1MaxSpace,RAW2MaxSpace,RAW3MaxSpace,RAW4MaxSpace,needSaveFile,needJuliaSend,Logfolder)
# println(cfgObj)
# cfgObj

max_stars_repo_name: wangzhongtian/fiberAlarmAnalyseTool

max_issues_repo_name: wangzhongtian/fiberAlarmAnalyseTool

max_forks_repo_name: wangzhongtian/fiberAlarmAnalyseTool


text:
function simulate_sim_random_walk(model::SingleIndexModel, initial_price::Float64,
    market_return::Array{Float64,1}; N::Int64 = 100, Δt = 1.0)

    # get the number of steps -
    number_of_steps = length(market_return)

    # initialize some storage -
    simulated_price_array = Array{Float64,2}(undef, number_of_steps, N)
    simulated_price_array[1,1:N] .= initial_price

    # get stuff from the model -
    α = model.α
    β = model.β
    r = model.r
    ϵ = model.ϵ

    # main loop -
    for path_index ∈ 1:N

        # what is the random pertrubation?
        ϵ_array = rand(ϵ, number_of_steps)

        # simulation loop -
        for simulation_step_index ∈ 2:number_of_steps

            # what is the return for the market at this time index?
            r_market = market_return[simulation_step_index - 1]

            # ok, compute the return in this time period -
            μ = α + β*(r_market - r) + ϵ_array[simulation_step_index] - r

            # now that we have the return μ, let's calc the new price -
            old_price = simulated_price_array[simulation_step_index - 1, path_index]
            new_price = old_price*(1+μ*Δt)
            simulated_price_array[simulation_step_index, path_index] = new_price
        end
    end

    # return -
    return simulated_price_array
end

function simulate_insample_portfolio_allocation(tickers::Array{String,1},
    return_dictionary::Dict{String,DataFrame}, ω::Array{Float64,1}, initial_budget::Float64,
    start::Date, stop::Date; multiplier::Float64 = 1.0)

    # initialize -
    df = extract_data_block_for_date_range(return_dictionary[tickers[1]], start, stop);
    number_of_timesteps = nrow(df)
    number_of_tickers = length(tickers)
    wealth_array = Array{Float64,2}(undef, number_of_timesteps, number_of_tickers)

    # initialize the wealth array -
    for ticker_index ∈ 1:number_of_tickers
        wealth_array[1,ticker_index] = ω[ticker_index]*initial_budget
    end

    for (index, ticker) ∈ enumerate(tickers)

        # grab -
        df = extract_data_block_for_date_range(return_dictionary[ticker], start, stop);

        for time_index ∈ 2:number_of_timesteps

            # get the return -
            r = multiplier*df[time_index, :μ];

            # compute the wealth -
            old_wealth = wealth_array[time_index-1,index]
            new_wealth = (1+r)*old_wealth
            wealth_array[time_index,index] = new_wealth
        end
    end

    # return the wealth_array -
    return wealth_array
end

function simulate_random_walk_model_trajectory(model::Distribution, initial_price::Float64,
    number_of_steps::Int64; number_of_sample_paths = 1, number_of_strata = 1)

    # initialize -
    number_of_steps = number_of_steps + 1
    sample_return_data = Array{Array{Float64,1},1}(undef, number_of_steps)

	# Let's use stratefied sampling to generate the return samples -
    for time_step_index ∈ 1:number_of_steps

		tmp_vector = Array{Float64,1}()

        # sample the strata ...
        for strata_index ∈ 1:number_of_strata

            # compute a number_of_sample_paths draws from tis strata?
            for _ ∈ 1:number_of_sample_paths

                # compute V -
                V₁ = (strata_index - 1)/number_of_strata + rand()/number_of_strata
                V₂ = (strata_index - 1)/number_of_strata + (1-rand())/number_of_strata

                # compute the quantile for this V -
                q₁ = quantile(model, V₁)
                q₂ = quantile(model, V₂)

                # grab this value -
                push!(tmp_vector, q₁)
                push!(tmp_vector, q₂)
            end
        end

        sample_return_data[time_step_index] = tmp_vector
    end

    # how many samples do we need?
    #sample_return_array = rand(model, number_of_steps, number_of_sample_paths)
    sample_return_array = transpose(hcat(sample_return_data...))

    # ok - so how many samples did we generate?
    (_, number_of_total_samples) = size(sample_return_array)
    price_array = Array{Float64,2}(undef, number_of_steps, number_of_total_samples)
    price_array[1, 1:number_of_total_samples] .= initial_price

    # compute the price -
    for sample_path_index = 1:number_of_total_samples
        for step_index = 2:number_of_steps
            price_array[step_index, sample_path_index] = price_array[step_index-1, sample_path_index] + sample_return_array[step_index, sample_path_index]
        end
    end

    # return -
    return price_array
end

# short cut methods -
(model::SingleIndexModel)(initial_price::Float64, market_return::Array{Float64,1}; Δt = 1.0) =
    simulate_sim_random_walk(model, initial_price, market_return; Δt = Δt);

# short cut method RWMC method -
(model::Distribution)(initial_price::Float64, number_of_steps::Int64;
    number_of_sample_paths = 1, number_of_strata = 1) =
    simulate_random_walk_model_trajectory(model, initial_price, number_of_steps;
        number_of_sample_paths = number_of_sample_paths, number_of_strata = number_of_strata)

max_stars_repo_name: Pooksoft/Serenity

max_issues_repo_name: Alba-Intelligence/Serenity

max_forks_repo_name: Alba-Intelligence/Serenity


text:
using Base.Test
using QuantumOptics

ωc = 1.2
ωa = 0.9
g = 1.0
γ = 0.5
κ = 1.1

T = Float64[0.,1.]

fockbasis = FockBasis(10)
spinbasis = SpinBasis(1//2)
basis = tensor(spinbasis, fockbasis)

sx = sigmax(spinbasis)
sy = sigmay(spinbasis)
sz = sigmaz(spinbasis)
sp = sigmap(spinbasis)
sm = sigmam(spinbasis)

Ha = embed(basis, 1, 0.5*ωa*sz)
Hc = embed(basis, 2, ωc*number(fockbasis))
Hint = sm ⊗ create(fockbasis) + sp ⊗ destroy(fockbasis)
H = Ha + Hc + Hint
Hdense = full(H)

Ja = embed(basis, 1, sqrt(γ)*sm)
Ja2 = embed(basis, 1, sqrt(0.5*γ)*sp)
Jc = embed(basis, 2, sqrt(κ)*destroy(fockbasis))
J = [Ja, Jc]
Jdense = map(full, J)

Ψ₀ = spinup(spinbasis) ⊗ fockstate(fockbasis, 2)
ρ₀ = Ψ₀⊗dagger(Ψ₀)

tout, ρt = timeevolution.master([0,100], ρ₀, Hdense, Jdense; reltol=1e-7)

ρss = steadystate.master(Hdense, Jdense; eps=1e-4)
@test tracedistance(ρss, ρt[end]) < 1e-3

ρss = steadystate.eigenvector(Hdense, Jdense)
@test tracedistance(ρss, ρt[end]) < 1e-6

ρss = steadystate.eigenvector(H, J)
@test tracedistance(ρss, ρt[end]) < 1e-3


text:
module Restful

using HttpServer
using URIParser

include("constants.jl")
include("RestMethod.jl")
include("Resource.jl")
include("Mixin.jl")
include("route.jl")
include("macros.jl")

include("utils/hooks.jl")
include("utils/mixins.jl")
include("utils/staticserver.jl")
include("utils/template.jl")

end # module Restful


text:
module TestExpressionsFusion

import Test: @testset, @test
import Volcanito: fuse_conjunction

@testset "fuse_conjunction" begin
    test_cases = (
        (
            input = (:(x > y || z == 0), :(y < z)),
            output = :((x > y || z == 0) && y < z),
        ),
        (
            input = (:(x), :(y)),
            output = :(x && y),
        ),
        (
            input = (:(x > y), :(y < z)),
            output = :(x > y && y < z),
        ),
        (
            input = (:(a > 0), :(b < 1)),
            output = :(a > 0 && b < 1),
        ),
    )
    for case in test_cases
        @test fuse_conjunction(case.input) == case.output
    end
end

end


text:
include(joinpath(@__DIR__, "testcommon.jl"))

# Wrap the functions defined in C++
module CppInheritance

using CxxWrap
@wrapmodule(Main.libinheritance, :define_types_module)

function __init__()
  @initcxx
end

export A, B, C, D, message, create_abstract, shared_ptr_message, shared_b, shared_c, shared_d, weak_ptr_message_a, weak_ptr_message_b, dynamic_message_c, take_ref

end

module VirtualSolver
  using CxxWrap
  @wrapmodule(Main.libinheritance, :define_vsolver_module)

  function __init__()
    @initcxx
  end
end

using .CppInheritance

@testset "$(basename(@__FILE__)[1:end-3])" begin

b = B()
c = C()
global d = D()

@test message(b) == "B"
@test message(c) == "C"
@test message(d) == "D"

@test take_ref(b) == "B"
@test take_ref(c) == "C"
@test take_ref(d) == "D"

# factory function returning an abstract type A
@test message(create_abstract()) == "B"

@test dynamic_message_c(c) == "C"

# shared ptr variants
b_ptr = shared_b()
c_ptr = shared_c()
d_ptr = shared_d()

@test shared_ptr_message(b_ptr) == "B"
@test shared_ptr_message(c_ptr) == "C"
@test shared_ptr_message(d_ptr) == "D"

@test message(b_ptr) == "B"
@test message(c_ptr) == "C"
@test message(d_ptr) == "D"

@test weak_ptr_message_b(b_ptr) == "B"
@test weak_ptr_message_a(b_ptr) == "B"
@test weak_ptr_message_a(c_ptr) == "C"
@test weak_ptr_message_a(d_ptr) == "D"

a = VirtualSolver.E()
VirtualSolver.solve(a)

b = VirtualSolver.F(@safe_cfunction(x -> 2x, Float64, (Float64,)))
VirtualSolver.solve(b)

end


text:
using GitForge
using GitForge: AStr, Endpoint, ForgeAPINotImplemented, ForgeAPIError
const GF = GitForge

using HTTP: HTTP
using JSON2: JSON2
using Test: @test, @testset, @test_throws

function capture(f::Function)
    so = stdout
    r, w = redirect_stdout()
    out = @async read(r, String)
    result = f()
    redirect_stdout(so)
    close(w)
    return result, fetch(out)
end

struct TestPostProcessor <: GF.PostProcessor end
GF.postprocess(::TestPostProcessor, ::HTTP.Response, ::Type) = :foo

struct TestForge <: GF.Forge end
GF.base_url(::TestForge) = "https://httpbin.org"
GF.request_headers(::TestForge, ::Function) = ["Foo" => "Bar"]
GF.request_query(::TestForge, ::Function) = Dict("foo" => "bar")
GF.request_kwargs(::TestForge, ::Function) = Dict(:verbose => 1)
GF.postprocessor(::TestForge, ::Function) = TestPostProcessor()
GF.endpoint(::TestForge, ::typeof(get_user)) = GF.Endpoint(:GET, "/get")
GF.into(::TestForge, ::typeof(get_user)) = Symbol

@testset "GitForge.jl" begin
    f = TestForge()
    (val, resp), out = capture(() -> get_user(f))

    @testset "Basics" begin
        @test val === :foo
        @test resp isa HTTP.Response
    end

    @testset "Request options" begin
        body = JSON2.read(IOBuffer(resp.body))
        @test startswith(get(body, :url, ""), "https://httpbin.org")
        @test get(get(body, :headers, Dict()), :Foo, "") == "Bar"
        @test get(get(body, :args, Dict()), :foo, "") == "bar"
        @test occursin("/get", get(body, :url, ""))
        @test !isempty(out)
    end

    @testset "Per-call request options" begin
        (_, resp), out = capture() do
            get_user(
                f;
                query=Dict("a" => "b"),
                headers=["A" => "B"],
                request_opts=(; verbose=0),
            )
        end

        @test isempty(out)

        body = JSON2.read(IOBuffer(resp.body))
        @test haskey(body.headers, :Foo)
        @test get(body.headers, :A, "") == "B"
        @test haskey(get(body, :args, Dict()), :foo)
        @test get(get(body, :args, Dict()), :a, "") == "b"
    end
end

# test whether apis are conformant

function test_api_func(api, (func, args...), auth)
    exclude(api, func, args...) && return
    try
        endpoint = GF.endpoint(api, func, args...)
        headers = Dict(GF.request_headers(api, func))
        @test headers["Content-Type"] isa AStr
        @test match(r"^Julia v.* \(GitForge v.*\)", headers["User-Agent"]) !== nothing
        if auth
            @test headers["Authorization"] isa AStr
        end
        @test GF.has_rate_limits(api, func) isa Bool
        @test !GF.has_rate_limits(api, func) || GF.on_rate_limit(api, func) isa GF.OnRateLimit
        @test endpoint isa Endpoint
        if hasmethod(GF.into, typeof.((api, func)))
            @test GF.into(api, func) isa Type
        else
            @test typeof(GF.postprocessor(api, func)) in [GF.DoNothing, GF.DoSomething, JSON]
        end
        true
    catch err
        !(err isa ForgeAPINotImplemented) && rethrow(err)
        err isa ForgeAPINotImplemented && !err.noted && println(sprint(io-> showerror(io, err)))
        err isa ForgeAPINotImplemented && !err.noted && @test true skip = true
        false
    end
end

mock_id(api) = mock_id(user_id_type(api))
mock_id(::Type{Integer}) = 1
mock_id(::Type{String}) = "fred"

function test_api(api; auth = false)
    unimplemented = 0
    @test GF.base_url(api) isa AStr
    for func in [
        get_user, (get_user, "fred"), (get_user, mock_id(api)),
        get_users,
        update_user, (update_user, mock_id(api)),
        create_user, (delete_user, mock_id(api)),
        get_user_repos, (get_user_repos, "fred"), (get_user_repos, mock_id(api)),
        (get_repo, "owner/repo"), (get_repo, "owner", "repo"), (get_repo, "owner", "subgroup", "repo"),
        create_repo, (create_repo, "org"), (create_repo, mock_id(api)), (create_repo, "namespace", "repo",
        (get_branch, "owner", "repo", "branch"),
        (get_branches, "owner", "repo"),
        (delete_branch, "owner", "repo", "branch"),
        (get_file_contents, "owner", "repo", "path"), (get_file_contents, mock_id(api), "path"),
        (get_pull_request, "owner", "repo", mock_id(api)), (get_pull_request, mock_id(api), mock_id(api)),
        (get_pull_requests, "owner", "repo"), (get_pull_requests, mock_id(api)),
        (create_pull_request, "owner", "repo"), (create_pull_request, mock_id(api)),
        (update_pull_request, "owner", "repo", 1), (update_pull_request, mock_id(api), 1),
        (subscribe_to_pull_request, mock_id(api), mock_id(api)),
        (unsubscribe_from_pull_request, mock_id(api), mock_id(api)),
        (get_commit, "owner", "repo", "ref"), (get_commit, mock_id(api), "ref"),
        (is_collaborator, "owner", "repo", "user"), (is_collaborator, "owner", "repo", mock_id(api)),
        (is_member, "org", "user"), (is_member, "org", mock_id(api)),
        groups,
        (get_tags, "owner", "repo"), (get_tags, mock_id(api)),
        (list_pull_request_comments, "owner", "repo", mock_id(api)), (list_pull_request_comments, mock_id(api), mock_id(api)),
        (get_pull_request_comment, "owner", "repo", mock_id(api)), (get_pull_request_comment, mock_id(api), mock_id(api), mock_id(api)),
        (create_pull_request_comment, "owner", "repo", mock_id(api)), (create_pull_request_comment, mock_id(api), mock_id(api)),
        (update_pull_request_comment, "owner", "repo", mock_id(api)), (update_pull_request_comment, mock_id(api), mock_id(api), mock_id(api)),
        (delete_pull_request_comment, "owner", "repo", mock_id(api)), (delete_pull_request_comment, mock_id(api), mock_id(api), mock_id(api)),
        (list_pipeline_schedules, mock_id(api)), (list_pipeline_schedules, "owner", "repo"),
        ]
        tuple = func isa Function ? (func,) : func
        impl = test_api_func(api, tuple, auth)
        !impl && (unimplemented += 1)
    end
    unimplemented > 0 && println("    $unimplemented unimplmented method$(unimplemented == 1 ? "" : "s")")
end

using GitForge.Bitbucket
using GitForge.GitHub
using GitForge.GitHub: Token, GitHubAPI
using GitForge.GitLab
using GitForge.GitLab: GitLabAPI
using UUIDs: UUID

const GH_TOKEN = Token("secret token")
const GL_TOKEN = OAuth2Token("secret token")
const USER = "a user"
const WORKSPACE = "a workspace"

user_id_type(::GitHubAPI) = Integer
user_id_type(::GitLabAPI) = Integer
user_id_type(::BitbucketAPI) = UUID

mock_id(::Type{UUID}) = UUID(0)

"""
    exclude(api, func, args...)

exclude an api function from testing
"""
exclude(_api, _func, _args...) = false
#exclude(::BitbucketAPI, ::typeof(get_file_contents), _args...) = true

@testset "Bitbucket.jl" begin
    api = BitbucketAPI(; token = GH_TOKEN, workspace = WORKSPACE)
    @test api.token == GH_TOKEN
    @test api.url == Bitbucket.DEFAULT_URL
    @test api.workspace == WORKSPACE
    test_api(api, auth = true)
    @test_throws ForgeAPIError GF.endpoint(api, get_repo, "owner")
end

@testset "GitHub.jl" begin
    api = GitHubAPI(; token = GH_TOKEN)
    @test api.token == GH_TOKEN
    @test api.url == GitHub.DEFAULT_URL
    test_api(api, auth = true)
end

@testset "GitLab.jl" begin
    api = GitLabAPI(; token = GL_TOKEN)
    @test api.token == GL_TOKEN
    @test api.url == GitLab.DEFAULT_URL
    test_api(api, auth = true)
end


text:
# This file contains a type for
# interacting with the generator

"""
Module containing the type to interact with the time series generator.
"""
module Gen

# Export
export GenStruct, gen!

# Load extras
using ..Extras

# Input data
using ..Input

# Arrays with non-standard indexing
using OffsetArrays

# Precisions and formats of numbers
using ..Prec

# Random numbers
using Random

"""
    GenStruct()

Instantiate this type to interact with the time series generator.

# Data
- `N::`[`IT`](@ref)`=0`: sample size;
- `Δt::`[`RT`](@ref)`=0`: sample step;
- `q::`[`RT`](@ref)`=0`: significance level;
- `α::`[`RT`](@ref)`=0`: parameter α of a linear trend;
- `β::`[`RT`](@ref)`=0`: parameter β of a linear trend;
- `r::`[`IT`](@ref)`=0`: number of harmonics;
- `A::Vector{`[`RT`](@ref)`}=[]`: amplitudes array;
- `ν::Vector{`[`RT`](@ref)`}=[]`: frequencies array;
- `ϕ::Vector{`[`RT`](@ref)`}=[]`: phase shifts array;
- `γ::`[`RT`](@ref)`=0`: «signal-to-noise» ratio.

# Methods
- [`read!`](@ref)`(file::AbstractString)`: read generator parameters from a file;
- [`example`](@ref)`(file::AbstractString)`: generate an example of a file containing
  the generator parameters;
- [`gen!`](@ref)`()`: generate time series;
- [`reset!`](@ref)`()`: reset an instance to default values.

# Note
Data can be also read when calling an instance like so:
```jldoctest; output = false
using Scats
s = Scats.API()
file, _ = mktemp()
s.Gen.example(file)
s.Gen(file)

# output

```
"""
mutable struct GenStruct

    # Data
    N::IT         # Sample size
    Δt::RT        # Sample step
    q::RT         # Significance level
    α::RT         # Parameter α of a linear trend
    β::RT         # Parameter β of a linear trend
    r::IT         # Number of harmonics
    A::Vector{RT} # Amplitudes array
    ν::Vector{RT} # Frequencies array
    ϕ::Vector{RT} # Phase shifts array
    γ::RT         # «Signal-to-noise» ratio

    # Methods
    read!::Function   # Read generator parameters from a file
    example::Function # Generate an example of a file containing the generator parameters
    gen!::Function    # Generate time series
    reset!::Function  # Reset an instance to default values

    # Construct an object of this type
    function GenStruct()
        this = new(0, 0, 0, 0, 0, 0, [], [], [], 0)
        this.read! = function (file::AbstractString) read!(this, file) end
        this.example = example
        this.gen! = function (Gen::GenStruct, Input::InputStruct) gen!(this, Input) end
        this.reset! = function () reset!(this) end
        return this
    end

    # Read the parameters when calling an instance
    function (Gen::GenStruct)(file::AbstractString)
        Gen.read!(file)
    end

end

# Sources
include("gen_exceptions.jl") # Exceptions
include("gen_read.jl")       # Read generator parameters from a file
include("gen_example.jl")    # Generate an example of a file containing
                             # the generator parameters
include("gen_gen.jl")        # Generate time series
include("gen_reset.jl")      # Reset an instance to default values

end


text:
#  Method
#  1. Argument reduction: Reduce x to an r so that |r| <= 0.5*log10(2). Given x,
#     find r and integer k such that
#
#                x = k*log10(2) + r,  |r| <= 0.5*log10(2).
#
# 2. Approximate exp10(r) by a polynomial on the interval [-0.5*log10(2), 0.5*log10(2)]:
#
#           exp10(x) = 1.0 + polynomial(x),
#
#    sup norm relative error within the interval of the polynomial approximations:
#    Float64 : [2.7245504724394698952e-18; 2.7245529895753476720e-18]
#    Float32 : [9.6026471477842205871e-10; 9.6026560194009888672e-10]
#
# 3. Scale back: exp10(x) = 2^k * exp10(r)

# log2(10)
const LOG2_10 = 3.321928094887362347870319429489390175864831393024580612054756395815934776608624
# log10(2)
const LOG10_2 = 3.010299956639811952137388947244930267681898814621085413104274611271081892744238e-01
# log(10)
const LN10 = 2.302585092994045684017991454684364207601101488628772976033327900967572609677367

# log10(2) into upper and lower bits
LOG10_2U(::Type{Float64}) = 3.01025390625000000000e-1
LOG10_2U(::Type{Float32}) = 3.00781250000000000000f-1

LOG10_2L(::Type{Float64}) = 4.60503898119521373889e-6
LOG10_2L(::Type{Float32}) = 2.48745663981195213739f-4

# max and min arguments
MAX_EXP10(::Type{Float64}) = 3.08254715559916743851e2 # log 2^1023*(2-2^-52)
MAX_EXP10(::Type{Float32}) = 38.531839419103626f0     # log 2^127 *(2-2^-23)

# one less than the min exponent since we can sqeeze a bit more from the exp10 function
MIN_EXP10(::Type{Float64}) = -3.23607245338779784854769e2 # log10 2^-1075
MIN_EXP10(::Type{Float32}) = -45.15449934959718f0         # log10 2^-150

@inline exp10_kernel(x::Float64) =
    @horner(x, 1.0,
    2.30258509299404590109361379290930926799774169921875,
    2.6509490552391992146397114993305876851081848144531,
    2.03467859229323178027470930828712880611419677734375,
    1.17125514891212478829629617393948137760162353515625,
    0.53938292928868392106522833273629657924175262451172,
    0.20699584873167015119932443667494226247072219848633,
    6.8089348259156870502017966373387025669217109680176e-2,
    1.9597690535095281527677713029333972372114658355713e-2,
    5.015553121397981796436571499953060992993414402008e-3,
    1.15474960721768829356725927226534622604958713054657e-3,
    1.55440426715227567738830671828509366605430841445923e-4,
    3.8731032432074128681303432086835414338565897196531e-5,
    2.3804466459036747669197886523306806338950991630554e-3,
    9.3881392238209649520573607528461934634833596646786e-5,
    -2.64330486232183387018679354696359951049089431762695e-2)

@inline exp10_kernel(x::Float32) =
    @horner(x, 1.0f0,
    2.302585124969482421875f0,
    2.650949001312255859375f0,
    2.0346698760986328125f0,
    1.17125606536865234375f0,
    0.5400512218475341796875f0,
    0.20749187469482421875f0,
    5.2789829671382904052734375f-2)

@eval exp10_small_thres(::Type{Float64}) = $(2.0^-29)
@eval exp10_small_thres(::Type{Float32}) = $(2.0f0^-14)

"""
    exp10(x)

Compute ``10^x``.

# Examples
```jldoctest
julia> exp10(2)
100.0

julia> exp10(0.2)
1.5848931924611136
```
"""
function exp10(x::T) where T<:Union{Float32,Float64}
    xa = reinterpret(Unsigned, x) & ~sign_mask(T)
    xsb = signbit(x)

    # filter out non-finite arguments
    if xa > reinterpret(Unsigned, MAX_EXP10(T))
        if xa >= exponent_mask(T)
            xa & significand_mask(T) != 0 && return T(NaN)
            return xsb ? T(0.0) : T(Inf) # exp10(+-Inf)
        end
        x > MAX_EXP10(T) && return T(Inf)
        x < MIN_EXP10(T) && return T(0.0)
    end
    # compute approximation
    if xa > reinterpret(Unsigned, T(0.5)*T(LOG10_2)) # |x| > 0.5 log10(2).
        # argument reduction
        if xa < reinterpret(Unsigned, T(1.5)*T(LOG10_2)) # |x| <= 1.5 log10(2)
            if xsb
                k = -1
                r = LOG10_2U(T) + x
                r = LOG10_2L(T) + r
            else
                k = 1
                r = x - LOG10_2U(T)
                r = r - LOG10_2L(T)
            end
        else
            n = round(T(LOG2_10)*x)
            k = unsafe_trunc(Int,n)
            r = muladd(n, -LOG10_2U(T), x)
            r = muladd(n, -LOG10_2L(T), r)
        end
        # compute approximation on reduced argument
        y = exp10_kernel(r)
        # scale back
        if k > -significand_bits(T)
            # multiply by 2.0 first to prevent overflow, extending the range
            k == exponent_max(T) && return y * T(2.0) * T(2.0)^(exponent_max(T) - 1)
            twopk = reinterpret(T, rem(exponent_bias(T) + k, fpinttype(T)) << significand_bits(T))
            return y*twopk
        else
            # add significand_bits(T) + 1 to lift the range outside the subnormals
            twopk = reinterpret(T, rem(exponent_bias(T) + significand_bits(T) + 1 + k, fpinttype(T)) << significand_bits(T))
            return y * twopk * T(2.0)^(-significand_bits(T) - 1)
        end
    elseif xa < reinterpret(Unsigned, exp10_small_thres(T))  # |x| < exp10_small_thres
        # Taylor approximation for small values: exp10(x) ≈ 1.0 + log(10)*x
        return muladd(x, T(LN10), T(1.0))
    else
        # primary range with k = 0, so compute approximation directly
        return exp10_kernel(x)
    end
end

max_stars_repo_name: nextjournal/julia-lang

max_issues_repo_name: nextjournal/julia-lang

max_forks_repo_name: nextjournal/julia-lang


text:
#=
These functions are the reason KeyedArray's keys are Ref not Tuple when ndims==1.
They should always change A.keys[] to have the right length,
but if the type doesn't fit, then return something not === to original
=#

"""
    push!(A::KeyedArray, val)

This adds `val` to the end of `pareent(A)`, and attempts to extend `axiskeys(A,1)` by one.
"""
function Base.push!(A::KeyedArray, val)
    data = push!(parent(A), val)
    old_keys = axiskeys(A,1)
    new_keys = extend_one!!(old_keys)
    if new_keys === old_keys  # then it was mutated, best case
        return A
    elseif typeof(new_keys) <: typeof(old_keys)  # will fit into same container
        setindex!(A.keys, new_keys)
        return A
    else  # we can't fix A, but can return correct thing
        return KeyedArray(data, (new_keys,))
    end
end

function Base.pop!(A::KeyedArray)
    val = pop!(parent(A))
    old_keys = axiskeys(A,1)
    new_keys = shorten_one!!(old_keys)
    if new_keys === old_keys
    elseif typeof(new_keys) <: typeof(old_keys)
        setindex!(A.keys, new_keys)
    else
        error("failed to shorten keys of array")
    end
    val
end

extend_one!!(r::Base.OneTo) = Base.OneTo(last(r)+1)
extend_one!!(r::UnitRange{Int}) = UnitRange(r.start, r.stop + 1)
extend_one!!(r::StepRange{<:Any,Int}) = StepRange(r.start, r.step, r.stop + r.step)
extend_one!!(r::Vector{<:Number}) = push!(r, length(r)+1)
extend_one!!(r::AbstractVector) = vcat(r, length(r)+1)

shorten_one!!(r::Base.OneTo) = Base.OneTo(last(r)-1)
shorten_one!!(r::Vector) = pop!(r)
shorten_one!!(r::AbstractVector) = r[1:end-1]

function Base.append!(A::KeyedArray, B)
    data = append!(parent(A), B)
    old_keys = axiskeys(A,1)
    new_keys = extend_by!!(old_keys, length(B))
    if new_keys === old_keys
        return A
    elseif typeof(new_keys) <: typeof(old_keys)
        setindex!(A.keys, new_keys)
        return A
    else
        return KeyedArray(data, (new_keys,))
    end
end
function Base.append!(A::KeyedArray, B::KeyedArray)
    data = append!(parent(A), parent(B))
    old_keys = axiskeys(A,1)
    new_keys = append!!(old_keys, axiskeys(B,1))
    if new_keys === old_keys
        return A
    elseif typeof(new_keys) <: typeof(old_keys)
        setindex!(A.keys, new_keys)
        return A
    else
        return KeyedArray(data, (new_keys,))
    end
end

extend_by!!(r::Base.OneTo, n::Int) = Base.OneTo(last(r)+n)
extend_by!!(r::UnitRange{Int}, n::Int) = UnitRange(r.start, r.stop + n)
extend_by!!(r::StepRange{<:Any,Int}, n::Int) = StepRange(r.start, r.step, r.stop + n * r.step)
extend_by!!(r::Vector{<:Number}, n::Int) = append!(r, length(r)+1 : length(r)+n+1)
extend_by!!(r::AbstractVector, n::Int) = vcat(r, length(r)+1 : length(r)+n+1)

append!!(r::Vector, s::AbstractVector) = append!(r,s)
append!!(r::AbstractVector, s::AbstractVector) = (extend_by!!(r, length(s)); vcat(r,s))

"""
    push!(A::KeyedArray; key = val)
    push!(A::KeyedArray, key => val)

This pushes `val` into `A.data`, and pushes `key` to `axiskeys(A,1)`.
Both of these must be legal operations, e.g. `A = wrapdims([1], ["a"]); push!(A, b=2)`.
"""
Base.push!(A::KeyedArray; kw...) = push!(A, map(Pair, keys(kw), values(values(kw)))...)

function Base.push!(A::KeyedArray, pairs::Pair...)
    axiskeys(A,1) isa AbstractRange && error("can't use push!(A, key => val) when axiskeys(A,1) isa AbstractRange")
    T = eltype(axiskeys(A,1))
    push!(axiskeys(A,1),  map(p -> T(first(p)), pairs)...)
    push!(parent(A), map(last, pairs)...)
    A
end


text:
# specialized methods
"""
    dotdot(::Vec, ::SymmetricFourthOrderTensor, ::Vec)

Computes a special dot product between two vectors and a symmetric fourth order tensor
such that ``a_k C_{ikjl} b_l``.
"""
@generated function dotdot(v1::Vec{dim}, S::SymmetricTensor{4, dim}, v2::Vec{dim}) where {dim}
    idx(i,j,k,l) = compute_index(SymmetricTensor{4, dim}, i, j, k, l)
    exps = Expr(:tuple)
    for j in 1:dim, i in 1:dim
        ex1, ex2 = Expr[],Expr[]
        for l in 1:dim, k in 1:dim
            push!(ex1, :(get_data(v1)[$k] * get_data(S)[$(idx(i,k,j,l))]))
            push!(ex2, :(get_data(v2)[$l]))
        end
        push!(exps.args, reducer(ex1, ex2, true))
    end
    return quote
        $(Expr(:meta, :inline))
        @inbounds return Tensor{2, dim}($exps)
    end
end


text:
# Generative model w/ Glow architecture from Kingma & Dhariwal (2018)
# Author: Philipp Witte, pwitte3@gatech.edu
# Date: January 2020

using InvertibleNetworks, LinearAlgebra, Flux

import Flux.Optimise.update!

# Define network
nx = 64     # must be multiple of 2
ny = 64
n_in = 3
n_hidden = 32
batchsize = 10
L = 2   # number of scales
K = 2   # number of flow steps per scale

# Input
X = rand(Float32, nx, ny, n_in, batchsize) #|> gpu

# Glow network
G = NetworkGlow(n_in, n_hidden, L, K) #|>gpu

# Objective function
function loss(X)
    Y, logdet = G.forward(X)
    f = .5f0/batchsize*norm(Y)^2 - logdet
    ΔX, X_ = G.backward(1f0./batchsize*Y, Y)
    return f
end

# Evaluate loss
f = loss(X)

# Update weights
opt = Flux.ADAM()
Params = get_params(G)
for p in Params
    update!(opt, p.data, p.grad)
end
clear_grad!(G)