text:
__precompile__()

"""
This is a package I use to handle numerical-model parameters,
thus the name. However, it should be useful otherwise too.
It has two main features:
- keyword type constructors with default values, and
- unpacking and packing of composite types and dicts.

The macro `@with_kw` which decorates a type definition to
allow default values and a keyword constructor:

```
julia> using Parameters

julia> @with_kw struct A
           a::Int = 6
           b::Float64 = -1.1
           c::UInt8
       end

julia> A(c=4)
A
  a: 6
  b: -1.1
  c: 4
```

Unpacking is done with `@unpack` (`@pack!` is similar):
```
struct B
    a
    b
    c
end
@unpack a, c = B(4,5,6)
# is equivalent to
BB = B(4,5,6)
a = BB.a
c = BB.c
```
"""
module Parameters
import Base: @__doc__
import OrderedCollections: OrderedDict
using UnPack: @unpack, @pack!

export @with_kw, @with_kw_noshow, type2dict, reconstruct, @unpack, @pack!, @pack

## Parser helpers
#################

# To iterate over code blocks dropping the line-number bits:
struct Lines
    block::Expr
end
start(lns::Lines) = 1
function next(lns::Lines, nr)
    for i=nr:length(lns.block.args)
        if lns.block.args[i] isa LineNumberNode
            continue
        end
        if ( lns.block.args[i] isa Symbol
             || lns.block.args[i] isa String # doc-string
             || !(lns.block.args[i].head==:line))
            return lns.block.args[i], i+1
        end
    end
    return -1
end
function done(lns::Lines, nr)
    if next(lns::Lines, nr)==-1
        true
    else
        false
    end
end
function Base.iterate(lns::Lines, nr=start(lns))
    nr = next(lns, nr)
    return nr == -1 ? nothing : nr
end

# This is not O(1) but hey...
function Base.setindex!(lns::Lines, val, ind)
    ii = 1
    for i=1:length(lns.block.args)
        if lns.block.args[i] isa LineNumberNode
            continue
        end
        if lns.block.args[i] isa Symbol || !(lns.block.args[i].head==:line)
            if ind==ii
                lns.block.args[i] = val
                return nothing
            end
            ii +=1
        end
    end
    throw(BoundsError("Attempted to set line $ind of $(ii-1) length code-block"))
end

# Transforms :(a::b) -> :a
decolon2(a::Expr) = (@assert a.head==:(::);  a.args[1])
decolon2(a::Symbol) = a

# Keep the ::T of the args if T ∈ typparas,
# leave symbols as is, drop field-doc-strings.
function keep_only_typparas(args, typparas)
    args = copy(args)
    tokeep = Int[]
    typparas_ = map(stripsubtypes, typparas)
    for i=1:length(args)
        isa(args[i], String) && continue # do not keep field doc-strings
        push!(tokeep, i)
        isa(args[i], Symbol) && continue
        # keep the typepara if ∈ typparas
        @assert args[i].head==:(::)
        T = args[i].args[2]
        if !(symbol_in(typparas_, T))
            args[i] = decolon2(args[i])
        end
    end
    args[tokeep]
end

# check whether a symbol is contained in an expression
symbol_in(s::Symbol, ex::Symbol) = s==ex
symbol_in(s::Symbol, ex) = false
function symbol_in(s::Symbol, ex::Expr)
    for a in ex.args
        symbol_in(s,a) && return true
    end
    return false
end
symbol_in(s::Symbol, ex::Vector) = any(map(e->symbol_in(s, e), ex))
symbol_in(s::Vector, ex) = any(map(ss->symbol_in(ss, ex), s))

# Returns the name of the type as Symbol
function typename(typedef::Expr)
    if typedef.args[2] isa Symbol
        return typedef.args[2]
    elseif typedef.args[2].args[1] isa Symbol
        return typedef.args[2].args[1]
    elseif typedef.args[2].args[1].args[1] isa Symbol
        return typedef.args[2].args[1].args[1]
    else
        error("Could not parse type-head from: $typedef")
    end
end

# Transforms:  Expr(:<:, :A, :B) -> :A
stripsubtypes(s::Symbol) = s
function stripsubtypes(e::Expr)
    e.args[1]
end
stripsubtypes(vec::Vector) = [stripsubtypes(v) for v in vec]

function check_inner_constructor(l)
    if l.args[1].head==:where
        fnhead = l.args[1].args[1]
    else
        fnhead = l.args[1]
    end
    if length(fnhead.args)==1
        error("No inner constructors with zero positional arguments allowed!")
    elseif (length(fnhead.args)==2 #1<length(fnhead.args)<=3
            && fnhead.args[2] isa Expr
            && fnhead.args[2].head==:parameters)
        error("No inner constructors with zero positional arguments plus keyword arguments allowed!")
    end
    nothing
end

## Exported helper functions
#####################
"""
Transforms a type-instance into a dictionary.

```
julia> type T
           a
           b
       end

julia> type2dict(T(4,5))
Dict{Symbol,Any} with 2 entries:
  :a => 4
  :b => 5
```

Note that this uses `getproperty`.
"""
function type2dict(dt)
    di = Dict{Symbol,Any}()
    for n in propertynames(dt)
        di[n] = getproperty(dt, n)
    end
    di
end

"""
Make a new instance of a type with the same values as the input type
except for the fields given in the AbstractDict second argument or as
keywords.  Works for types, Dicts, and NamedTuples.

Note: this is not very performant.  Check Setfield.jl for a faster &
nicer implementation.

```jldoctest
julia> using Parameters

julia> struct A
           a
           b
       end

julia> a = A(3,4)
A(3, 4)

julia> b = reconstruct(a, b=99)
A(3, 99)
```
"""
function reconstruct(pp::T, di) where T
    if pp isa AbstractDict
        pp = deepcopy(pp)
        for (k,v) in di
            !haskey(pp, k) && error("Field $k not in type $T")
            pp[k] = v
        end
        return pp
    else
        di = !isa(di, AbstractDict) ? Dict(di) : copy(di)
        ns = fieldnames(T)
        args = []
        for (i,n) in enumerate(ns)
            push!(args, pop!(di, n, getfield(pp, n)))
        end
        length(di)!=0 && error("Fields $(keys(di)) not in type $T")
        return pp isa NamedTuple ? T(Tuple(args)) : T(args...)
    end
end
reconstruct(pp; kws...) = reconstruct(pp, kws)

###########################
# Keyword constructors with @with_kw
##########################

# A type with fields (r,a) in variable aa becomes
# quote
#     r = aa.r
#     a = aa.a
# end
_unpack(binding, fields) = Expr(:block, [:($f = $binding.$f) for f in fields]...)
# Pack fields back into binding using reconstruct:
function _pack_mutable(binding, fields)
    e = Expr(:block, [:($binding.$f = $f) for f in fields]...)
    push!(e.args, binding)
    e
end
function _pack_new(T, fields)
    Expr(:call, T, fields...)
end

const macro_hidden_nargs = length(:(@m).args) - 1 # ==1 on Julia 0.6, ==2 on Julia 0.7

"""
This function is called by the `@with_kw` macro and does the syntax
transformation from:

```julia
@with_kw struct MM{R}
    r::R = 1000.
    a::R
end
```

into

```julia
struct MM{R}
    r::R
    a::R
    MM{R}(r,a) where {R} = new(r,a)
    MM{R}(;r=1000., a=error("no default for a")) where {R} = MM{R}(r,a) # inner kw, type-paras are required when calling
end
MM(r::R,a::R) where {R} = MM{R}(r,a) # default outer positional constructor
MM(;r=1000,a=error("no default for a")) =  MM(r,a) # outer kw, so no type-paras are needed when calling
MM(m::MM; kws...) = reconstruct(mm,kws)
MM(m::MM, di::Union{AbstractDict, Tuple{Symbol,Any}}) = reconstruct(mm, di)
macro unpack_MM(varname)
    esc(quote
    r = varname.r
    a = varname.a
    end)
end
macro pack_MM(varname)
    esc(quote
    varname = $Parameters.reconstruct(varname,r=r,a=a)
    end)
end
```
"""
function with_kw(typedef, mod::Module, withshow=true)
    if typedef.head==:tuple # named-tuple
        withshow==false && error("`@with_kw_noshow` not supported for named tuples")
        return with_kw_nt(typedef, mod)
    elseif typedef.head != :struct
        error("""Only works on type-defs or named tuples.
              Make sure to have a space after `@with_kw`, e.g. `@with_kw (a=1,)
              Also, make sure to use a trailing comma for single-field NamedTuples.
              """)
    end
    err1str = "Field \'"
    err2str = "\' has no default, supply it with keyword."

    inner_constructors = Any[]

    # parse a few things
    tn = typename(typedef) # the name of the type
    ismutable = typedef.args[1]
    # Returns M{...} (removes any supertypes)
    if typedef.args[2] isa Symbol
        typparas = Any[]
    elseif typedef.args[2].head==:<:
        if typedef.args[2].args[1] isa Symbol
            typparas = Any[]
        else
            typparas = typedef.args[2].args[1].args[2:end]
        end
    else
        typparas = typedef.args[2].args[2:end]
    end

    # error on types without fields
    lns = Lines(typedef.args[3])
    if done(lns, start(lns))
        error("@with_kw only supported for types which have at least one field.")
    end
    # default type @deftype
    l, i = next(lns, start(lns))
    if l isa Expr && l.head == :macrocall && l.args[1] == Symbol("@deftype")
        has_deftyp = true
        if length(l.args) != (2 + macro_hidden_nargs)
            error("Malformed `@deftype` line $l")
        end
        deftyp = l.args[2 + macro_hidden_nargs]
        if done(lns, i)
            error("@with_kw only supported for types which have at least one field.")
        end
    else
        has_deftyp = false
    end

    # Expand all macros (except @assert) in body now (only works at
    # top-level)
    # See issue https://github.com/mauro3/Parameters.jl/issues/21
    lns2 = Any[] # need new lines as expanded macros may have many lines
    for (i,l) in enumerate(lns) # loop over body of typedef
        if i==1 && has_deftyp
            push!(lns2, l)
            continue
        end
        if l isa Symbol || l isa String
            push!(lns2, l)
            continue
        end
        if l.head==:macrocall && l.args[1]!=Symbol("@assert")
            tmp = macroexpand(mod, l)
            if tmp.head==:block
                llns = Lines(tmp)
                for ll in llns
                    push!(lns2, ll)
                end
            else
                push!(lns2,tmp)
            end
        else
            push!(lns2, l)
        end
    end
    lns = lns2

    # the vars for the unpack macro
    unpack_vars = Any[]
    # the type def
    fielddefs = quote end # holds r::R etc
    fielddefs.args = Any[] # in julia 0.5 this is [:( # /home/mauro/.julia/v0.5/Parameters/src/Parameters.jl, line 228:)]
    kws = OrderedDict{Any, Any}()
    # assertions in the body
    asserts = Any[]
    for (i,l) in enumerate(lns) # loop over body of typedef
        if i==1 && has_deftyp
            # ignore @deftype line
            continue
        end
        if l isa Symbol  # no default value and no type annotation
            if has_deftyp
                push!(fielddefs.args, :($l::$deftyp))
            else
                push!(fielddefs.args, l)
            end
            sym = l
            syms = string(sym)
            kws[sym] = :(error($err1str * $syms * $err2str))
            # unwrap-macro
            push!(unpack_vars, sym)
        elseif l isa String # doc-string
            push!(fielddefs.args, l)
        elseif l.head==:(=)  # default value and with or without type annotation
            if l.args[1] isa Expr && (l.args[1].head==:call || # inner constructor
                                        l.args[1].head==:where && l.args[1].args[1].head==:call) # inner constructor with `where`
                check_inner_constructor(l)
                push!(inner_constructors, l)
            else
                fld = l.args[1]
                if fld isa Symbol && has_deftyp # no type annotation
                    fld = :($fld::$deftyp)
                end
                # add field doc-strings
                docstring = string("Default: ", l.args[2])
                if i > 1 && lns[i-1] isa String
                    # if the last line was a docstring, append the default
                    fielddefs.args[end] *= " " * docstring
                else
                    # otherwise add a new line
                    push!(fielddefs.args, docstring)
                end
                push!(fielddefs.args, fld)
                kws[decolon2(fld)] = l.args[2]
                # unwrap-macro
                push!(unpack_vars, decolon2(fld))
            end
        elseif l.head==:macrocall  && l.args[1]==Symbol("@assert")
            # store all asserts
            push!(asserts, l)
        elseif l.head==:function # inner constructor
            check_inner_constructor(l)
            push!(inner_constructors, l)
        elseif l.head==:block
            error("No nested begin-end allowed in type defintion")
        else # no default value but with type annotation
            push!(fielddefs.args, l)
            sym = decolon2(l.args[1])
            syms = string(sym)
            kws[sym] = :(error($err1str *$syms * $err2str))
            # unwrap-macro
            push!(unpack_vars, l.args[1])
        end
    end
    # The type definition without inner constructors:
    typ = Expr(:struct, typedef.args[1:2]..., copy(fielddefs))

    # Inner keyword constructor.  Note that this calls the positional
    # constructor under the hood and not `new`.  That way a user can
    # provide a special positional constructor (say enforcing
    # invariants) which also gets used with the keywords.
    args = Any[]
    kwargs = Expr(:parameters)
    for (k,w) in kws
        push!(args, k)
        push!(kwargs.args, Expr(:kw,k,w))
    end
    if length(typparas)>0
        tps = stripsubtypes(typparas)
        innerc = :( $tn{$(tps...)}($kwargs) where {$(tps...)} = $tn{$(tps...)}($(args...)))
    else
        innerc = :($tn($kwargs) = $tn($(args...)) )
    end
    push!(typ.args[3].args, innerc)

    # Inner positional constructor: only make it if no inner
    # constructors are user-defined.  If one or several are defined,
    # assume that one has the standard positional signature.
    if length(inner_constructors)==0
        if length(typparas)>0
            tps = stripsubtypes(typparas)
            innerc2 = :( $tn{$(tps...)}($(args...)) where {$(tps...)} = new{$(tps...)}($(args...)) )
        else
            innerc2 = :($tn($(args...)) = new($(args...)))
        end
        prepend!(innerc2.args[2].args, asserts)
        push!(typ.args[3].args, innerc2)
    else
        if length(asserts)>0
            error("Assertions are only allowed in type-definitions with no inner constructors.")
        end
        append!(typ.args[3].args, inner_constructors)
    end

    # Outer positional constructor which does not need explicit
    # type-parameters when called.  Only make this constructor if
    #  (1) type parameters are used at all
    #  (2) all type parameters are used in the fields (otherwise get a
    #      "method is not callable" warning!)
    #       See also https://github.com/JuliaLang/julia/issues/17186
    if typparas!=Any[] # condition (1)
        # fields definitions stripped of ::Int etc., only keep ::T if T∈typparas :
        fielddef_strip_contT = keep_only_typparas(fielddefs.args, typparas)
        outer_positional = :(  $tn($(fielddef_strip_contT...)) where {$(typparas...)}
                             = $tn{$(stripsubtypes(typparas)...)}($(args...)))
        # Check condition (2)
        checks = true
        for tp in stripsubtypes(typparas)
            checks = checks && symbol_in(tp, fielddefs.args)
        end
        if !checks
            outer_positional = :()
        end
    else
        outer_positional = :()
    end

    # Outer keyword constructor, useful to infer the type parameter
    # automatically.  This calls the outer positional constructor.
    # only create if type parameters are used.
    if typparas==Any[]
        outer_kw=:()
    else
        outer_kw = :($tn($kwargs) = $tn($(args...)) )
    end
    # NOTE: The reason to have both outer and inner keyword
    # constructors are to allow both calls:
    #   `MT4(r=4, a=5.0)` (outer kwarg-constructor) and
    #   `MT4{Float32, Int}(r=4, a=5.)` (inner kwarg constructor).

    ## outer copy constructor
    ###
    outer_copy = quote
        $tn(pp::$tn; kws... ) = $Parameters.reconstruct(pp, kws)
        # $tn(pp::$tn, di::Union(AbstractDict,Vararg{Tuple{Symbol,Any}}) ) = reconstruct(pp, di) # see issue https://github.com/JuliaLang/julia/issues/11537
        # $tn(pp::$tn, di::Union(AbstractDict, Tuple{Vararg{Tuple{Symbol, Any}}}) ) = reconstruct(pp, di) # see issue https://github.com/JuliaLang/julia/issues/11537
        $tn(pp::$tn, di::$Parameters.AbstractDict) = $Parameters.reconstruct(pp, di)
        $tn(pp::$tn, di::Vararg{Tuple{Symbol,Any}} ) = $Parameters.reconstruct(pp, di)
    end

    # (un)pack macro from https://groups.google.com/d/msg/julia-users/IQS2mT1ITwU/hDtlV7K1elsJ
    unpack_name = Symbol("unpack_"*string(tn))
    pack!_name = Symbol("pack_"*string(tn)*"!")
    pack_name = Symbol("pack_"*string(tn))
    showfn = if withshow
        :(function Base.show(io::IO, p::$tn)
              if get(io, :compact, false) || get(io, :typeinfo, nothing)==$tn
                Base.show_default(IOContext(io, :limit => true), p)
              else
                # just dumping seems to give ok output, in particular for big data-sets:
                dump(IOContext(io, :limit => true), p, maxdepth=1)
              end
          end)
    else
        :nothing
    end
    if ismutable
        pack_macros = quote
            macro $pack!_name(ex)
                esc($Parameters._pack_mutable(ex, $unpack_vars))
            end
            macro $pack_name()
                esc($Parameters._pack_new($tn, $unpack_vars))
            end
        end
    else
        pack_macros = quote
            macro $pack_name()
                esc($Parameters._pack_new($tn, $unpack_vars))
            end
        end
    end

    # Finish up
    quote
        Base.@__doc__ $typ
        $outer_positional
        $outer_kw
        $outer_copy
        $showfn
        macro $unpack_name(ex)
            esc($Parameters._unpack(ex, $unpack_vars))
        end
        $pack_macros
        $tn
    end
end

"""
Do the with-kw stuff for named tuples.
"""
function with_kw_nt(typedef, mod)
    kwargs = []
    args = []
    nt = []
    for a in typedef.args
        if a isa Expr
            a.head != :(=) && error("NameTuple fields need to be of form: `k=val`")
            sy = a.args[1]::Symbol
            va = a.args[2]
            push!(kwargs, Expr(:kw, sy, va))
            push!(args, sy)
            push!(nt, :($sy=$sy))
        elseif a isa Symbol  # no default value given
            sy = a
            push!(args, sy)
            push!(nt, :($sy=$sy))
            push!(kwargs, Expr(:kw, sy, :(error("Supply default value for $($(string(sy)))"))))
        else
            error("Cannot parse $(string(a))")
        end
    end
    NT = gensym(:NamedTuple_kw)
    nt = Expr(:tuple, nt...)
    quote
        $NT(; $(kwargs...)) =$nt
        $NT($(args...)) = $nt
        $NT
    end
end

"""
Macro which allows default values for field types and a few other features.

Basic usage:

```julia
@with_kw struct MM{R}
    r::R = 1000.
    a::Int = 4
end
```

For more details see manual.
"""
macro with_kw(typedef)
    return esc(with_kw(typedef, __module__, true))
end

macro with_kw(args...)
    error("""Only works on type-defs or named tuples.
          Did you try to construct a NamedTuple but omitted the space between the macro and the NamedTuple?
          Do `@with_kw (a=1, b=2)` and not `@with_kw(a=1, b=2)`.
          """)
end

"""
As `@with_kw` but does not define a `show` method to avoid annoying
redefinition warnings.

```julia
@with_kw_noshow struct MM{R}
    r::R = 1000.
    a::Int = 4
end
```

For more details see manual.
"""
macro with_kw_noshow(typedef)
    return esc(with_kw(typedef, __module__, false))
end

end # module


text:
using Printf
using KiteModels, KitePodModels, KiteUtils

# change this to KPS3 or KPS4
const Model = KPS4

if ! @isdefined kcu;  const kcu = KCU(se());   end
if ! @isdefined kps4; const kps4 = Model(kcu); end

# the following values can be changed to match your interest
dt = 0.05
STEPS = 600
PLOT = true
FRONT_VIEW = false
ZOOM = true
PRINT = false
STATISTIC = false
# end of user parameter section #

if PLOT
    using Plots
    include("plot2d.jl")
end

function simulate(integrator, steps, plot=false)
    start = integrator.p.iter
    for i in 1:steps
        if PRINT
            lift, drag = KiteModels.lift_drag(kps4)
            @printf "%.2f: " round(integrator.t, digits=2)
            println("lift, drag  [N]: $(round(lift, digits=2)), $(round(drag, digits=2))")
        end

        KiteModels.next_step!(kps4, integrator, dt=dt)

        if plot
            reltime = i*dt
            if mod(i, 5) == 0
                p = plot2d(kps4.pos, reltime; zoom=ZOOM, front=FRONT_VIEW, segments=se().segments)
                display(p)
            end
        end
    end
    (integrator.p.iter - start) / steps
end

integrator = KiteModels.init_sim!(kps4, stiffness_factor=0.04, prn=STATISTIC)

if PLOT
    av_steps = simulate(integrator, STEPS, true)
else
    println("\nStarting simulation...")
    simulate(integrator, 100)
    runtime = @elapsed av_steps = simulate(integrator, STEPS-100)
    println("\nTotal simulation time: $(round(runtime, digits=3)) s")
    speed = (STEPS-100) / runtime * dt
    println("Simulation speed: $(round(speed, digits=2)) times realtime.")
end
lift, drag = KiteModels.lift_drag(kps4)
println("lift, drag  [N]: $(round(lift, digits=2)), $(round(drag, digits=2))")
println("Average number of callbacks per time step: $av_steps")

ext: jl


text:
using SoapySDR

# Don't forget to add/import a device-specific plugin package!
# using xtrx_jll
# using SoapyLMS7_jll
# using SoapyRTLSDR_jll

for (didx, dev) in enumerate(Devices())
    @info("device", dev, idx=didx)

    @info("TX channels:")
    for (idx, tx_channel) in enumerate(dev.tx)
        display(tx_channel)
    end

    @info("RX channels:")
    for (idx, rx_channel) in enumerate(dev.rx)
        display(rx_channel)
    end
end


text:
function yaxis!(ec::EChart; formatter::Union{String, JSFunction, Nothing} = nothing, kwargs...)

	for (k, v) in kwargs
	   setfield!(ec.yAxis[1], k, v)
	end

	#Apply axis label format if present
	formatter != nothing ? ec.yAxis[1].axisLabel.formatter = formatter : nothing

	return ec

end

function xaxis!(ec::EChart; formatter::Union{String, JSFunction, Nothing} = nothing, kwargs...)

	for (k, v) in kwargs
	   setfield!(ec.xAxis[1], k, v)
	end

	#Apply axis label format if present
	formatter != nothing ? ec.xAxis[1].axisLabel.formatter = formatter : nothing

	return ec

end


text:
__precompile__(true)
module Foo
import Qux
this = "Foo2"
end

max_stars_repo_name: nreigl/julia

max_issues_repo_name: BenChung/Julia-BTI

max_forks_repo_name: BenChung/Julia-BTI


text:
# Random point in the Dim-dimensional unit hypercube
Base.rand(::Type{Point{Dim, T}}) where {Dim, T} = Point{Dim, T}(rand(SVector{Dim, T}))
function Base.rand(::Type{Point{Dim, T}}, num_points::Int64) where {Dim, T}
    return [ Point{Dim, T}(rand(SVector{Dim, T})) for i = 1:num_points ]
end

# Random line in the Dim-dimensional unit hypercube
function Base.rand(::Type{LineSegment{Dim, F}}) where {Dim, F}
    points = rand(Point{Dim, F}, 2)
    return LineSegment{Dim, F}(points[1], points[2])
end

# N random lines in the Dim-dimensional unit hypercube
function Base.rand(::Type{LineSegment{Dim, F}}, N::Int64) where {Dim, F}
    return [ rand(LineSegment{Dim, F}) for i ∈ 1:N ]
end

# Random quadratic segment in the Dim-dimensional unit hypercube
function Base.rand(::Type{QuadraticSegment{Dim, F}}) where {Dim, F}
    points = rand(Point{Dim, F}, 3)
    return QuadraticSegment(points[1], points[2], points[3])
end

# N random quadratic segments in the Dim-dimensional unit hypercube
function Base.rand(::Type{QuadraticSegment{Dim, F}}, N::Int64) where {Dim, F}
    return [ rand(QuadraticSegment{Dim, F}) for i ∈ 1:N ]
end

# A random AABox within the Dim-dimensional unit hypercube
# What does the distribution of AABoxs look like? Is this uniform?
function Base.rand(::Type{AABox{Dim, T}}) where {Dim, T}
    coord₁ = rand(T, Dim)
    coord₂ = rand(T, Dim)
    return AABox{Dim, T}(Point{Dim, T}(min.(coord₁, coord₂)),
                         Point{Dim, T}(max.(coord₁, coord₂)))
end

# N random AABoxs within the Dim-dimensional unit hypercube
function Base.rand(::Type{AABox{Dim, T}}, num_boxes::Int64) where {Dim, T}
    return [ rand(AABox{Dim, T}) for i ∈ 1:num_boxes ]
end

max_stars_repo_name: khurrumsaleem/MOCNeutronTransport

max_issues_repo_name: khurrumsaleem/MOCNeutronTransport

max_forks_repo_name: khurrumsaleem/MOCNeutronTransport


text:
const seen = Dict{Vector{Char}, Vector{Char}}()

function findnextterm(prevterm)
    counts = Dict{Char, Int}()
    reversed = Vector{Char}()
    for c in prevterm
        if !haskey(counts, c)
            counts[c] = 0
        end
        counts[c] += 1
    end
    for c in sort(collect(keys(counts)))
        if counts[c] > 0
            push!(reversed, c)
            if counts[c] == 10
                push!(reversed, '0'); push!(reversed, '1')
            else
                push!(reversed, Char(UInt8(counts[c]) + UInt8('0')))
            end
        end
    end
    reverse(reversed)
end

function findsequence(seedterm)
    term = seedterm
    sequence = Vector{Vector{Char}}()
    while !(term in sequence)
        push!(sequence, term)
        if !haskey(seen, term)
            nextterm = findnextterm(term)
            seen[term] = nextterm
        end
        term = seen[term]
    end
    return sequence
end

function selfseq(maxseed)
    maxseqlen = -1
    maxsequences = Vector{Pair{Int, Vector{Char}}}()
    for i in 1:maxseed
        seq = findsequence([s[1] for s in split(string(i), "")])
        seqlen = length(seq)
        if seqlen > maxseqlen
            maxsequences = [Pair(i, seq)]
            maxseqlen = seqlen
        elseif seqlen == maxseqlen
            push!(maxsequences, Pair(i, seq))
        end
    end
    println("The longest sequence length is $maxseqlen.")
    for p in maxsequences
        println("\n Seed: $(p[1])")
        for seq in p[2]
            println("     ", join(seq, ""))
        end
    end
end

selfseq(1000000)

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
"""
Finite difference calculation & machine learning of electric potential & electric field from charge
"""

using LinearAlgebra
using Plots
using Random
using Flux
Random.seed!(1)
using Main.EquivariantOperators # end users should omit Main.

# make grid
dims = 2
dx = 0.1
cell = dx * Matrix(I, dims, dims)
rmax = 1.0
grid = Grid(cell, rmax)

# make operators
rmin = 1e-9
rmax = sqrt(3)
ϕ = Op(r -> 1 / (4π * r), rmin, rmax, cell)
E = Op(r -> 1 / (4π * r^2), rmin, rmax, cell; l = 1)
▽ = Op(:▽, cell)

# put dipole charges
ρf = zeros(size(grid)..., 1)
put!(ρf, grid, [0.5 0.0; -0.5 0.0]', [1.0 -1]'')

# calculate fields
Ef = E(ρf)
ϕf = ϕ(ρf)

# test
rvec = [0, 0]
@show get(ϕf, grid, rvec), [0.0]
@show get(Ef, grid, rvec), get(-▽(ϕ(ρf)), grid, rvec), [-2 / (4π * 0.5^2), 0]

p = []
push!(p, heatmap(ρf[:, :, 1]',title = "dipole charge"), )
push!(p, heatmap(ϕf[:, :, 1]', title = "dipole potential"),)
plot(p..., layout = length(p))
savefig("d1.svg");

# @unpack x, y=grid
x = grid.coords[:, :, 1]
y = grid.coords[:, :, 2]
s = 1e-2
u = s * Ef[:, :, 1]
v = s * Ef[:, :, 2]
x, y, u, v = vec.([x, y, u, v])
quiver(x, y, quiver = (u, v), c = :blue, title = "dipole electric field")
savefig("d2.svg");

##
# make neural operators
ϕ_ = Op(Radfunc(), rmin, rmax, cell)
E_ = Op(Radfunc(), rmin, rmax, cell; l = 1)

ps = Flux.params(ϕ_, E_)
function loss()
    remake!(E_)
    remake!(ϕ_)
    global E_f = E_(ρf)
    global ϕ_f = ϕ_(ρf)
    @show l = (nae(E_f, Ef) + nae(ϕ_f, ϕf)) / 2
    l
end

data = [()]
loss()
opt = ADAM(0.1)
Flux.@epochs 50 Flux.train!(loss, ps, data, opt)

## plot
p = []
push!(p, heatmap(ϕ_.kernel[:, :, 1], title = "learned potential kernel"))
r = 0:0.01:1
push!(
    p,
    plot(r, ϕ_.radfunc.(r), title = "learned potential kernel radial function"),
)
plot(p..., layout = length(p))
savefig("ml1.svg");

x = E_.grid.coords[:, :, 1]
y = E_.grid.coords[:, :, 2]
s = 1e-2
u = s * E_.kernel[:, :, 1]
v = s * E_.kernel[:, :, 2]
x, y, u, v = vec.([x, y, u, v])
quiver(x, y, quiver = (u, v), c = :blue, title = "learned E field kernel ")
savefig("ml2.svg");


text:
# Based on this guide: https://ase.tufts.edu/chemistry/lin/images/FortranMD_TeachersGuide.pdf

using Atomistic
using AtomsBase
using InteratomicPotentials
using NBodySimulator
using Plots
using Unitful
using UnitfulAtomic

N = 864
element = :Ar
box_size = 3.47786u"nm"
reference_temp = 94.4u"K"
thermostat_prob = 0.1 # this number was chosen arbitrarily
Δt = 1e-2u"ps"

initial_bodies = generate_bodies_in_cell_nodes(N, element, box_size, reference_temp)
initial_system = FlexibleSystem(initial_bodies, CubicPeriodicBoundaryConditions(austrip(box_size)))

eq_steps = 2000
eq_thermostat = AndersenThermostat(austrip(reference_temp), thermostat_prob / austrip(Δt))
eq_simulator = NBSimulator(Δt, eq_steps, thermostat = eq_thermostat)
potential = LennardJones(austrip(1.657e-21u"J"), austrip(0.34u"nm"))

eq_result = @time simulate(initial_system, eq_simulator, potential)

temp = @time plot_temperature(eq_result, 10)
energy = @time plot_energy(eq_result, 10)

prod_steps = 5000
prod_simulator = NBSimulator(Δt, prod_steps, t₀ = get_time_range(eq_result)[end])

prod_result = @time simulate(get_system(eq_result), prod_simulator, potential)

display(@time plot_temperature!(temp, prod_result, 10))
display(@time plot_energy!(energy, prod_result, 10))

# rdf = @time plot_rdf(prod_result, potential.σ, 0.05)
# display(rdf)
# savefig(rdf, "artifacts/argon_lj_rdf.svg")

;


text:
module TestInstantiate

using Test
using MathOptInterface
const MOI = MathOptInterface

struct DummyOptimizer <: MOI.AbstractOptimizer end
MOI.is_empty(::DummyOptimizer) = true

function _test_instantiate(T)
    function f()
        return MOI.Utilities.MockOptimizer(
            MOI.Utilities.UniversalFallback(MOI.Utilities.Model{T}()),
        )
    end
    optimizer_constructor =
        MOI.OptimizerWithAttributes(f, MOI.Silent() => true, "a" => 1, "b" => 2)
    optimizer = MOI.instantiate(optimizer_constructor)
    @test optimizer isa MOI.Utilities.MockOptimizer{
        MOI.Utilities.UniversalFallback{MOI.Utilities.Model{T}},
        Float64,
    }
    @test MOI.get(optimizer, MOI.Silent())
    optimizer = MOI.instantiate(optimizer_constructor, with_bridge_type = T)
    @test optimizer isa MOI.Bridges.LazyBridgeOptimizer{
        MOI.Utilities.MockOptimizer{
            MOI.Utilities.UniversalFallback{MOI.Utilities.Model{T}},
            Float64,
        },
    }
    @test MOI.get(optimizer, MOI.Silent())
    @test MOI.get(optimizer, MOI.RawOptimizerAttribute("a")) == 1
    @test MOI.get(optimizer, MOI.RawOptimizerAttribute("b")) == 2

    optimizer_constructor = MOI.OptimizerWithAttributes(DummyOptimizer, [])
    optimizer = MOI.instantiate(optimizer_constructor)
    @test optimizer isa DummyOptimizer
    optimizer = MOI.instantiate(optimizer_constructor, with_bridge_type = T)
    @test optimizer isa MOI.Bridges.LazyBridgeOptimizer{
        MOI.Utilities.CachingOptimizer{
            DummyOptimizer,
            MOI.Utilities.UniversalFallback{MOI.Utilities.Model{T}},
        },
    }

    err = ErrorException(
        "The provided `optimizer_constructor` returned a non-empty optimizer.",
    )
    function g()
        model = f()
        MOI.add_variable(model)
        return model
    end
    @test_throws err MOI.instantiate(g)
    @test_throws err MOI.instantiate(g, with_bridge_type = T)
    optimizer_constructor = MOI.OptimizerWithAttributes(g)
    @test_throws err MOI.instantiate(optimizer_constructor)
    @test_throws err MOI.instantiate(
        optimizer_constructor,
        with_bridge_type = T,
    )

    err = ErrorException(MOI._INSTANTIATE_NOT_CALLABLE_MESSAGE)
    @test_throws err MOI.OptimizerWithAttributes(1)
    @test_throws err MOI.OptimizerWithAttributes(1, MOI.Silent() => true)
    @test_throws err MOI.instantiate(1)
    @test_throws err MOI.instantiate(1, with_bridge_type = T)

    err = ErrorException(
        "The provided `optimizer_constructor` returned an object of type " *
        "$Int. Expected a MathOptInterface.AbstractOptimizer.",
    )
    h() = 1
    @test_throws err MOI.instantiate(h)
    @test_throws err MOI.instantiate(h, with_bridge_type = T)
    optimizer_constructor = MOI.OptimizerWithAttributes(h)
    @test_throws err MOI.instantiate(optimizer_constructor)
    @test_throws err MOI.instantiate(
        optimizer_constructor,
        with_bridge_type = T,
    )
end

function test_instantiate_Float64()
    _test_instantiate(Float64)
    return
end

function test_instantiate_Int()
    _test_instantiate(Int)
    return
end

function runtests()
    for name in names(@__MODULE__; all = true)
        if startswith("$name", "test_")
            @testset "$(name)" begin
                getfield(@__MODULE__, name)()
            end
        end
    end
end

end

TestInstantiate.runtests()


text:
# Trying to append dataframe and save it as matfile
using DataFrames
using Dates
using MAT
using FileTools, DataFrameTools
using Revise
using JuliaToolsForMagTIP
dir = raw"D:\GeoMag (main)\GeoMag_1"
allpaths = filelistall(r"\.mat",dir);
dtstr = gettag.("dt", allpaths)
dt = Date.(dtstr, dateformat"yyyymmdd")
stn = gettag.("stn", allpaths)
typ = gettag.("type", allpaths)

df = DataFrame(:filename => basename.(allpaths), :date => dt, :fullpath=>allpaths, :station=>stn, :color => 1, :type => typ)
filter!(:type => x -> occursin(r"(full|Full)", x), df)

sort!(df, :date)
unique!(df, :filename)

fulldt = df.date[1]:Day(1):df.date[end] |> collect
gd = groupby(df, :station)

dfs = []
for dfi in gd
    fulldf = DataFrame(:date=> fulldt);
    push!(dfs, outerjoin(dfi, fulldf, on=:date))
end

df2 = vcat(dfs...)
dropmissing!(df2, :station)
if !isequal(sort(df2,[:date,:station]), df)
    error("They should be identical")
end

## ^^^^^^^^^Above code is the same as those in data_report.jl^^^^^^^^^^
function MAT.matopen(fname::AbstractString, mode::AbstractString; compress::Bool = false)
    mode == "r"  ? matopen(fname, true , false, false, false, false, false)    :
    mode == "r+" ? matopen(fname, true , true , false, false, false, compress) :
    mode == "w"  ? matopen(fname, false, true , true , true , false, compress) :
    mode == "w+" ? matopen(fname, true , true , true , true , false, compress) :
    mode == "a"  ? matopen(fname, false, true , true , false, true, compress)  :
    mode == "a+" ? matopen(fname, true , true , true , false, true, compress)  :
    throw(ArgumentError("invalid open mode: $mode"))
end

df3 = filter(:date => x -> (Date(2017,10,1) > x > Date(2010,9,30)) ,df2)
gd_stn = groupby(df3, [:station])
targetfolder = raw"D:\GeoMag (main)\GeoMag4HMM"
extrema(df3.date) # (Date("2010-10-01"), Date("2017-09-30"))
for df_i in gd_stn
    stnm = df_i.station[1]
    fname = "D_$(stnm)_ti1.mat"
    fpath = joinpath(targetfolder, fname);
    if isfile(fpath)
        continue
    end
    sort!(df_i, :date)
    try
        target = matopen(fpath, "w"; compress=true)
        df000 = DataFrame();
        for path in df_i.fullpath
            file = matopen(path)
            M = read(file, "M")
            close(file)
            dfm = DataFrame(M[:,1:7],
                    [:year, :month,:day,:hour,:minute,:second,:data])

            filter!(AsTable(:) => (nt -> isgeomagdatavalid(nt)) , dfm)

            select!(dfm,
            [:year, :month,:day,:hour,:minute,:second] => ByRow((yyyy,mm,dd,hh,MM,ss) -> toordinal(yyyy,mm,dd,hh,MM,ss)) => :ordinaldate, :data)

            append!(df000, dfm)

        end
        write(target, "M", Matrix(df000))
        close(target)
    catch e
        rm(fpath);
        @warn "Error occurred while loading $path, writing to $fpath. $fpath has been deleted."
        rethrow(e)
    end
end


text:
using Primes

factortodict(d, n) = (d[minimum(collect(keys(factor(n))))] = n)

# Numbers are from from the Perl 6 example.
numbers = [64921987050997300559,  70251412046988563035,  71774104902986066597,
           83448083465633593921,  84209429893632345702,  87001033462961102237,
           87762379890959854011,  89538854889623608177,  98421229882942378967,
           259826672618677756753, 262872058330672763871, 267440136898665274575,
           278352769033314050117, 281398154745309057242, 292057004737291582187]

mins = Dict()

Base.@sync(
    Threads.@threads for n in numbers
        factortodict(mins, n)
    end
)

answer = maximum(keys(mins))
println("The number that has the largest minimum prime factor is $(mins[answer]), with a smallest factor of $answer")


text:
type ThresholdModel <: DischargeModel
    M::Int
    k::Int
    f::Function
    β::Matrix{Float64}
    H::Matrix{Float64}
    Q::Vector{Float64}
end

ThresholdModel(M,N,f,k=1) = ThresholdModel(M,k,f,zeros(M,k),zeros(M,N),zeros(N))

function clusters(M::ThresholdModel,H::Matrix{Float64})
    cs = zeros(Int,size(H,2))
    for i in 1:size(H,2)
        cs[i] = M.f(H[:,i])
    end
    cs
end

function estimate_by_clusters(H,Q,cs,M,k)
    β = zeros(M,k)
    for i in 1:k
        Hi = H[:,cs.==i]
        Qi = Q[cs.==i]
        β[:,i] = (Hi*Hi')\(Hi*Qi)
    end
    β
end

function StatsBase.fit(M::ThresholdModel,H::Matrix{Float64},Q::Vector{Float64})
    cs = clusters(M,H)
    β = estimate_by_clusters(H,Q,cs,M.M,M.k)
    ThresholdModel(M.M,M.k,M.f,β,H,Q)
end

function StatsBase.fit!(M::ThresholdModel,H::Matrix{Float64},Q::Vector{Float64})
    cs = clusters(M,H)
    β = estimate_by_clusters(H,Q,cs,M.M,M.k)
    M.β = β
    M.H = H
    M.Q = Q
    M
end

StatsBase.coef(M::ThresholdModel) = M.β
StatsBase.nobs(M::ThresholdModel) = length(M.Q)
StatsBase.model_response(M::ThresholdModel) = M.Q

function StatsBase.predict(M::ThresholdModel,H)
    cs = clusters(M,H)
    Q = zeros(size(H,2))
    for i in 1:length(Q)
        Q[i] = (H[:,i]'M.β[:,cs[i]])[1]
    end
    Q
end

StatsBase.predict(M::ThresholdModel) = predict(M,M.H)

StatsBase.residuals(M::ThresholdModel) = M.Q-predict(M)
StatsBase.residuals(M::ThresholdModel,H,Q) = Q-predict(M,H)


text:
using Gadfly
using RDatasets
using DataFrames

iris = data("datasets", "iris")

hstack(
    [plot(subset(iris, :(Species .== $(species))),
         x="Sepal.Length", y="Sepal.Width", Geom.point)
     for species in ["setosa", "versicolor", "virginica"]])


text:
using Test, LegibleLambdas

@testset "Basics" begin
    f = @λ(x -> x + 1)
    @test f isa Function
    @test repr(f) == "(x -> x + 1)"
    @test f(2) == 3

    g = @λ((x,y) -> x^2 + y^2)
    @test g isa Function
    @test repr(g) == "((x, y) -> x ^ 2 + y ^ 2)"
    @test g(-15.0, 10.0) == 325.0

    h = @λ((x, y, z) -> occursin(x, y*z))
    @test h isa Function
    @test repr(h) == "((x, y, z) -> occursin(x, y * z))"
    @test h("hi", "abh", "ing")
end

@testset "Closures" begin
    D(f, ϵ=1e-10) = @λ(x -> (f(x+ϵ)-f(x))/ϵ)

    if VERSION < v"1.1.0"
        @test repr(D(sin)) == "(x -> (f(x + ϵ) - f(x)) / ϵ)"
    else
        @test repr(D(sin)) == "(x -> ((sin)(x + 1.0e-10) - (sin)(x)) / 1.0e-10)"
    end

    @test D(sin)(π) == -1.000000082740371
end

@testset "Unicode" begin
    @test repr(@λ((η) -> η^2 + 1)) == "(η -> η ^ 2 + 1)"
end

@testset "Varargs" begin
    f = @λ((x=1) -> x^2 + 1)
    @test repr(f) == "(x = 1 -> x ^ 2 + 1)"

    g = @λ((x, y=1) -> x^2 + y)
    @test repr(g) == "((x, y = 1) -> x ^ 2 + y)"

    h = @λ((x, y=1, z=2) -> x^2 + y/z)
    @test repr(h) == "((x, y = 1, z = 2) -> x ^ 2 + y / z)"
end


text:
using LinearAlgebra

struct ValueOne; end
ValueOne()

# Compute X <- a X + b I.
function matfun_axpby!(X,a,b,Y::UniformScaling)
    m,n=size(X)
    if ~(a isa ValueOne)
        rmul!(X,a)
    end
    @inbounds for i=1:n
        X[i,i]+=(b isa ValueOne) ? 1 : b
    end
end

@inline function exp_ps_m11(A)
    T=promote_type(eltype(A),Float64)
    A_copy=similar(A,T); A_copy .= A;
    return exp_ps_m11!(A_copy)
end

@inline function exp_ps_m11!(A)
    T=promote_type(eltype(A),Float64) # Make it work for many 'bigger' types (matrices and scalars)
    # max_memslots=9
    n=size(A,1)
    # The first slots are precomputed nodes [:A]
    memslots2 = similar(A,T)
    memslots3 = similar(A,T)
    memslots4 = similar(A,T)
    memslots5 = similar(A,T)
    memslots6 = similar(A,T)
    memslots7 = similar(A,T)
    memslots8 = similar(A,T)
    memslots9 = similar(A,T)
    # Assign precomputed nodes memslots
    memslots1=A # overwrite A
    # Uniform scaling is exploited.
    # No matrix I explicitly allocated.
    value_one=ValueOne()
    # Computation order: A2 A3 A4 A5 A6 A7 P5 C4 P4 C3 P3 C2 P2 C1 P1 C0 P0
    # Computing A2 with operation: mult
    mul!(memslots2,memslots1,memslots1)
    # Computing A3 with operation: mult
    mul!(memslots3,memslots1,memslots2)
    # Computing A4 with operation: mult
    mul!(memslots4,memslots1,memslots3)
    # Computing A5 with operation: mult
    mul!(memslots5,memslots1,memslots4)
    # Computing A6 with operation: mult
    mul!(memslots6,memslots1,memslots5)
    # Computing A7 with operation: mult
    mul!(memslots7,memslots1,memslots6)
    # Computing P5 = x*I+x*A+x*A2+x*A3+x*A4+x*A5+x*A6+x*A7
    coeff1=9.67759295863189e-41
    coeff2=2.6882202662866363e-42
    coeff3=7.265460179153071e-44
    coeff4=1.911963205040282e-45
    coeff5=4.902469756513544e-47
    coeff6=1.2256174391283858e-48
    coeff7=2.9893108271424046e-50
    coeff8=7.117406731291439e-52
    memslots8 .= coeff2.*memslots1 .+ coeff3.*memslots2 .+ coeff4.*memslots3 .+ coeff5.*memslots4 .+ coeff6.*memslots5 .+ coeff7.*memslots6 .+ coeff8.*memslots7
    mul!(memslots8,true,I*coeff1,true,true)
    # Computing C4 with operation: mult
    mul!(memslots9,memslots8,memslots7)
    # Deallocating P5 in slot 8
    # Computing P4 = x*I+x*A+x*A2+x*A3+x*A4+x*A5+x*A6+x*C4
    coeff1=3.279889237069838e-30
    coeff2=1.1309962886447716e-31
    coeff3=3.7699876288159054e-33
    coeff4=1.216125041553518e-34
    coeff5=3.8003907548547434e-36
    coeff6=1.151633562077195e-37
    coeff7=3.387157535521162e-39
    coeff8=1.0
    # Smart lincomb recycle C4
    memslots9 .= coeff2.*memslots1 .+ coeff3.*memslots2 .+ coeff4.*memslots3 .+ coeff5.*memslots4 .+ coeff6.*memslots5 .+ coeff7.*memslots6 .+ coeff8.*memslots9
    mul!(memslots9,true,I*coeff1,true,true)
    # Computing C3 with operation: mult
    mul!(memslots8,memslots9,memslots7)
    # Deallocating P4 in slot 9
    # Computing P3 = x*I+x*A+x*A2+x*A3+x*A4+x*A5+x*A6+x*C3
    coeff1=1.9572941063391263e-20
    coeff2=8.896791392450574e-22
    coeff3=3.868170170630684e-23
    coeff4=1.6117375710961184e-24
    coeff5=6.446950284384474e-26
    coeff6=2.4795962632247976e-27
    coeff7=9.183689863795546e-29
    coeff8=1.0
    # Smart lincomb recycle C3
    memslots8 .= coeff2.*memslots1 .+ coeff3.*memslots2 .+ coeff4.*memslots3 .+ coeff5.*memslots4 .+ coeff6.*memslots5 .+ coeff7.*memslots6 .+ coeff8.*memslots8
    mul!(memslots8,true,I*coeff1,true,true)
    # Computing C2 with operation: mult
    mul!(memslots9,memslots8,memslots7)
    # Deallocating P3 in slot 8
    # Computing P2 = x*I+x*A+x*A2+x*A3+x*A4+x*A5+x*A6+x*C2
    coeff1=1.1470745597729725e-11
    coeff2=7.647163731819816e-13
    coeff3=4.779477332387385e-14
    coeff4=2.8114572543455206e-15
    coeff5=1.5619206968586225e-16
    coeff6=8.22063524662433e-18
    coeff7=4.110317623312165e-19
    coeff8=1.0
    # Smart lincomb recycle C2
    memslots9 .= coeff2.*memslots1 .+ coeff3.*memslots2 .+ coeff4.*memslots3 .+ coeff5.*memslots4 .+ coeff6.*memslots5 .+ coeff7.*memslots6 .+ coeff8.*memslots9
    mul!(memslots9,true,I*coeff1,true,true)
    # Computing C1 with operation: mult
    mul!(memslots8,memslots9,memslots7)
    # Deallocating P2 in slot 9
    # Computing P1 = x*I+x*A+x*A2+x*A3+x*A4+x*A5+x*A6+x*C1
    coeff1=0.0001984126984126984
    coeff2=2.48015873015873e-5
    coeff3=2.7557319223985893e-6
    coeff4=2.755731922398589e-7
    coeff5=2.505210838544172e-8
    coeff6=2.08767569878681e-9
    coeff7=1.6059043836821613e-10
    coeff8=1.0
    # Smart lincomb recycle C1
    memslots8 .= coeff2.*memslots1 .+ coeff3.*memslots2 .+ coeff4.*memslots3 .+ coeff5.*memslots4 .+ coeff6.*memslots5 .+ coeff7.*memslots6 .+ coeff8.*memslots8
    mul!(memslots8,true,I*coeff1,true,true)
    # Computing C0 with operation: mult
    mul!(memslots9,memslots8,memslots7)
    # Deallocating P1 in slot 8
    # Deallocating A7 in slot 7
    # Computing P0 = x*I+x*A+x*A2+x*A3+x*A4+x*A5+x*A6+x*C0
    coeff1=1.0
    coeff2=1.0
    coeff3=0.5
    coeff4=0.16666666666666666
    coeff5=0.041666666666666664
    coeff6=0.008333333333333333
    coeff7=0.001388888888888889
    coeff8=1.0
    # Smart lincomb recycle A
    memslots1 .= coeff2.*memslots1 .+ coeff3.*memslots2 .+ coeff4.*memslots3 .+ coeff5.*memslots4 .+ coeff6.*memslots5 .+ coeff7.*memslots6 .+ coeff8.*memslots9
    mul!(memslots1,true,I*coeff1,true,true)
    # Deallocating A2 in slot 2
    # Deallocating A3 in slot 3
    # Deallocating A4 in slot 4
    # Deallocating A5 in slot 5
    # Deallocating A6 in slot 6
    # Deallocating C0 in slot 9
    return memslots1 # Returning P0
end

max_stars_repo_name: matrixfunctions/GraphMatFunData

max_issues_repo_name: matrixfunctions/GraphMatFunData

max_forks_repo_name: matrixfunctions/GraphMatFunData


text:
B_data = [ -2.5   0.5        0.0   0.0   0.0   0.0   0.0   0.0   0.666667   0.0        0.0        0.0       0.0        0.0        0.0   0.0
1.0  -3.0        0.0   0.0   0.0   0.0   0.0   0.0   0.666667   0.0        0.0        0.0       0.0        0.0        0.0   0.0
0.0   0.0       -2.5   0.0   0.0   0.0   0.0   0.0   2.0        0.0        0.0        0.5       0.0        0.0        0.0   0.0
0.0   0.0        0.0  -2.5   0.0   0.0   0.0   0.0   0.0        2.0        0.0        0.0       0.5        0.0        0.0   0.0
0.0   0.0        0.0   0.0  -2.5   0.0   0.0   0.0   0.0        0.0        2.0        0.0       0.0        0.5        0.0   0.0
0.0   0.0        0.0   0.0   0.0  -2.5   0.0   0.0   2.0        0.0        0.0        0.5       0.0        0.0        0.0   0.0
0.0   0.0        0.0   0.0   0.0   0.0  -2.5   0.0   0.0        2.0        0.0        0.0       0.5        0.0        0.0   0.0
0.0   0.0        0.0   0.0   0.0   0.0   0.0  -2.5   0.0        0.0        2.0        0.0       0.0        0.5        0.0   0.0
0.0   0.0        1.0   0.0   0.0   0.0   0.0   0.0  -2.12821    0.435897  -0.961538   1.0       0.0        0.0        0.0   5.23077
0.0   0.0        0.0   1.0   0.0   0.0   0.0   0.0  -0.555556  -3.66667    3.38889    0.0       1.0        0.0        0.0  -9.33333
0.0   0.0        0.0   0.0   1.0   0.0   0.0   0.0   0.017094  -0.102564  -2.67735    0.0       0.0        1.0        0.0   6.10256
0.0   6.46154    0.0   0.0   0.0   1.0   0.0   0.0   2.0        0.0        0.0       -4.01603  -0.153846   0.025641   0.0   0.0
0.0  -4.0        0.0   0.0   0.0   0.0   1.0   0.0   0.0        2.0        0.0        5.08333  -5.5       -0.833333   0.0   0.0
0.0   0.538462   0.0   0.0   0.0   0.0   0.0   1.0   0.0        0.0        2.0       -1.44231   0.653846  -3.19231    0.0   0.0
0.0   0.0        0.0   0.0   0.0   0.0   0.0   0.0   0.0        0.0        0.0        0.0       0.0        0.0625    -2.5   2.0
0.0   0.0        0.0   0.0   0.0   0.0   0.0   0.0   0.0        0.0        0.0        0.0       0.0        0.125      1.0  -2.0]


text:
using Catlab
using Catlab.Doctrines
using Catlab.Graphics
using Catlab.WiringDiagrams
using Catlab.Programs
import Base.Multimedia: display
import Catlab.Graphics: to_graphviz, LeftToRight
import Base: (==), length, show
using Petri
using SemanticModels.CategoryTheory
import SemanticModels.CategoryTheory: undecorate, ⊔
using SemanticModels.PetriModels
using SemanticModels.PetriCospans
import SemanticModels.PetriCospans: otimes_ipm, compose_pushout
import SemanticModels: model

import Catlab.Doctrines:
    Ob,
    Hom,
    dom,
    codom,
    compose,
    ⋅,
    ∘,
    id,
    oplus,
    otimes,
    ⊗,
    ⊕,
    munit,
    mzero,
    braid,
    dagger,
    dunit,
    dcounit,
    mcopy,
    Δ,
    delete,
    ◊,
    mmerge,
    ∇,
    create,
    □,
    plus,
    zero,
    coplus,
    cozero,
    meet,
    top,
    join,
    bottom

export wd, draw, model, Disease

wd = to_wiring_diagram
draw(d::WiringDiagram) =
    to_graphviz(add_junctions(d), orientation = LeftToRight, labels = true)
draw(d::HomExpr) = draw(wd(d))

model(c::PetriCospan) = left(c.f).d[1].model

println("Making COVID Model")

@present Disease(FreeEpidemiology) begin
    S::Ob
    E::Ob
    I::Ob
    R::Ob
end


text:
using Millboard
using Test

Millboard.set_table_mode(:grid)
Base.eval(:(have_color = false))

Millboard.display_style[:prepend_newline] = true

board = [11 12 13; 21 22 23]
@test """

+---+----+----+----+
|   |  1 |  2 |  3 |
+===+====+====+====+
| 1 | 11 | 12 | 13 |
+---+----+----+----+
| 2 | 21 | 22 | 23 |
+---+----+----+----+""" == table(board) |> string

Millboard.display_style[:prepend_newline] = false

Base.eval(:(have_color = true))


text:
type Configuration{T <: Parameter} <: Parameter
    parameters::Dict{String, T}
    name::String
    value::Dict{String, Any}

    function Configuration(parameters::Dict{String, T}, name::String, values::Dict{String, Any})
        for key in keys(parameters)
            @inbounds values[key] = parameters[key].value
        end
        new(parameters, name, values)
    end
end

function Configuration{T <: Parameter}(parameters::Dict{String, T}, name::String)
    Configuration{T}(parameters, name, Dict{String, Any}())
end

function Configuration{T <: Parameter}(parameters::Array{T, 1}, name::String)
    params = Dict{String, T}()
    for parameter in parameters
        @inbounds params[parameter.name] = parameter
    end
    Configuration{T}(params, name, Dict{String, Any}())
end

function Configuration(name::String)
    params = Dict{String, Parameter}()
    Configuration{Parameter}(params, name, Dict{String, Any}())
end

function Base.convert{T <: Parameter}(::Type{Array{T}}, configuration::Configuration)
    parameter_array = T[]
    for key in collect(keys(configuration.parameters))
        push!(parameter_array, configuration[key])
    end
    parameter_array
end

function getindex(configuration::Configuration, index::String)
    configuration.parameters[index]
end

function setindex!{T <: Parameter}(configuration::Configuration, value::T, index::String)
    configuration.parameters[index] = value
end


text:
#------------------------------
#DEFINE STATIONARY SOURCE w/ OMNIDIRECTIONAL ANTENNA

struct LTISource
  position::Vector{Float64}
  transmission ::Function
end

# Method
function (𝚽::LTISource)(𝛏₀::Vector{Float64}, t₀::Float64)
   𝐩ₛ, p = 𝚽.position, 𝚽.transmission
   delay = distBetween(𝐩ₛ,𝛏₀)/lightSpeed
   return A(delay) * p(t₀-delay)
end

#-----------------------------------------------------------

struct Target
   coefficient::Vector{Float64}
   position::Vector{Vector{Float64}}
end

function(T::Target)(t₀::Float64)
αₖ = T.coefficient
𝛏ₖ = T.position

for i = 1:length(𝛏ₖ)
    r = αₖ[i].*𝚽(𝛏ₖ[i],t₀)  # transmission
end
end

max_stars_repo_path: src/Basic LTI Sources/sources.jl

max_issues_repo_path: src/Basic LTI Sources/sources.jl

max_forks_repo_path: src/Basic LTI Sources/sources.jl


text:
#!/usr/bin/env julia

using Luxor

using Test

using Random
Random.seed!(42)

function circle_as_path(fname)
    currentwidth = 1200
    currentheight = 1200 # pts
    Drawing(currentwidth, currentheight, fname)
    origin()
    background(setgray(0.8)...)
    tiles = Tiler(currentwidth, currentheight, 8, 8, margin=20)
    for (pos, n) in tiles
        randomhue()
        circlepath(pos, 70, :path)
        newsubpath()
        circlepath(pos, rand(5:65), :fill, reversepath=true)
    end
    @test finish() == true
end

fname = "circlepath.pdf"
circle_as_path(fname)
println("...finished test: output in $(fname)")


text:
import Base: size, ndims
import LinearAlgebra: diag

export ndoms,
       domainType,
       codomainType,
       is_linear,
       is_eye,
       is_null,
       is_diagonal,
       is_AcA_diagonal,
       is_AAc_diagonal,
       is_orthogonal,
       is_invertible,
       is_full_row_rank,
       is_full_column_rank,
       is_sliced,
       diag_AcA,
       diag_AAc,
       displacement,
       remove_displacement

"""
`domainType(A::AbstractOperator)`

Returns the type of the domain.

```julia
julia> domainType(DFT(10))
Float64

julia> domainType(hcat(Eye(Complex{Float64},(10,)),DFT(Complex{Float64},10)))
(Complex{Float64}, Complex{Float64})
```
"""
domainType

"""
`codomainType(A::AbstractOperator)`

Returns the type of the codomain.

```julia
julia> codomainType(DFT(10))
Complex{Float64}

julia> codomainType(vcat(Eye(Complex{Float64},(10,)),DFT(Complex{Float64},10)))
(Complex{Float64}, Complex{Float64})
```
"""
codomainType

"""
`size(A::AbstractOperator, [dom,])`

Returns the size of an `AbstractOperator`. Type `size(A,1)` for the size of the codomain and `size(A,2)` for the size of the codomain.
"""
size(L::AbstractOperator, i::Int) = size(L)[i]

"""
`ndims(A::AbstractOperator, [dom,])`

Returns a `Tuple` with the number of dimensions of the codomain and domain of an `AbstractOperator`.  Type `ndims(A,1)` for the number of dimensions of the codomain and `ndims(A,2)` for the number of dimensions of the codomain.

"""
ndims(L::AbstractOperator) = count_dims(size(L,1)), count_dims(size(L,2))
ndims(L::AbstractOperator, i::Int) = ndims(L)[i]

count_dims(dims::NTuple{N,Int}) where N = N
count_dims(dims::NTuple{N,Tuple}) where N = count_dims.(dims)

"""
`ndoms(L::AbstractOperator, [dom::Int]) -> (number of codomains, number of domains)`

Returns the number of codomains and domains  of a `AbstractOperator`. Optionally you can specify the codomain (with `dom = 1`) or the domain (with `dom = 2`)

```julia
julia > ndoms(DFT(10,10))
(1,1)

julia> ndoms(hcat(DFT(10,10),DFT(10,10)))
(1, 2)

julia> ndoms(hcat(DFT(10,10),DFT(10,10)),2)
2

julia> ndoms(DCAT(DFT(10,10),DFT(10,10)))
(2,2)
```
"""
ndoms(L::AbstractOperator) = length.(ndims(L))
ndoms(L::AbstractOperator, i::Int) = ndoms(L)[i]

diag_AcA(L::AbstractOperator) = error("is_AAc_diagonal($L) == false")
diag_AAc(L::AbstractOperator) = error("is_AcA_diagonal($L) == false")

"""
`is_linear(A::AbstractOperator)`

Test whether `A` is a `LinearOperator`

```julia
julia> is_linear(Eye(2))
true

julia> is_linear(Sigmoid(Float64,(2,)))
false
```

"""
is_linear(          L::LinearOperator  ) = true
is_linear(          L::AbstractOperator) = false

"""
`is_null(A::AbstractOperator)`

Test whether `A` is null.

```julia
julia> is_null(Zeros(Float64,(10,),(10,)))
true

julia> is_null(Eye(10))
false

```
"""
is_null(            L::AbstractOperator) = false

"""
`is_eye(A::AbstractOperator)`

Test whether `A` is an Identity operator

```julia
julia> is_eye(Eye(10))
true

julia> is_eye(Zeros(Float64,(10,),(10,)))
false

```
"""
is_eye(             L::AbstractOperator) = false

"""
`is_diagonal(A::AbstractOperator)`

Test whether `A` is diagonal.

```julia
julia> is_diagonal(Eye(10))
true

julia> is_diagonal(FiniteDiff((10,)))
false

```
"""
is_diagonal(        L::AbstractOperator) = false

"""
`is_AcA_diagonal(A::AbstractOperator)`

Test whether `A'*A` is diagonal.

```julia
julia> is_AcA_diagonal(Eye(10))
true

julia> is_AcA_diagonal(GetIndex((10,),1:3))
false

```
"""
is_AcA_diagonal(    L::AbstractOperator) = is_diagonal(L)

"""
`is_AAc_diagonal(A::AbstractOperator)`

Test whether `A*A'` is diagonal.

```julia
julia> is_AAc_diagonal(Eye(10))
true

julia> is_AAc_diagonal(GetIndex((10,),1:3))
false

```
"""
is_AAc_diagonal(    L::AbstractOperator) = is_diagonal(L)

"""
`is_orthogonal(A::AbstractOperator)`

Test whether `A` is orthogonal.

```julia
julia> is_orthogonal(DCT(10))
true

julia> is_orthogonal(MatrixOp(randn(3,4)))
false

```
"""
is_orthogonal(      L::AbstractOperator) = false

"""
`is_invertable(A::AbstractOperator)`

Test whether `A` is easily invertable.

```julia
julia> is_invertible(DFT(10))
true

julia> is_invertable(MatrixOp(randn(3,4)))
false

```
"""
is_invertible(      L::AbstractOperator) = false

"""
`is_full_row_rank(A::AbstractOperator)`

Test whether `A` is easily invertable.

```julia
julia> is_full_row_rank(MatrixOp(randn(3,4)))
true

julia> is_full_row_rank(MatrixOp(randn(4,3)))
false
```
"""
is_full_row_rank(   L::AbstractOperator) = false

"""
`is_full_row_rank(A::AbstractOperator)`

Test whether `A` is easily invertable.

```julia
julia> is_full_column_rank(MatrixOp(randn(4,3)))
true

julia> is_full_column_rank(MatrixOp(randn(3,4)))
false
```
"""
is_full_column_rank(L::AbstractOperator) = false

"""
`is_sliced(A::AbstractOperator)`

Test whether `A` is a sliced operator.

```julia
julia> is_sliced(GetIndex((10,), 1:5))
true

```
"""
is_sliced(L::AbstractOperator) = false

"""
`displacement(A::AbstractOperator)`

Returns the displacement of the operator.

```julia
julia> A = AffineAdd(Eye(4),[1.;2.;3.;4.])
I+d  ℝ^4 -> ℝ^4

julia> displacement(A)
4-element Array{Float64,1}:
 1.0
 2.0
 3.0
 4.0

```
"""
function displacement(S::AbstractOperator)
  D = domainType(S)
  if typeof(D) <: Tuple
    x = ArrayPartition(zeros.(D, size(S, 2))...)
  else
	  x = zeros(D, size(S, 2))
  end
  d = S*x
  if all(y -> y == d[1], d )
    return d[1]
  else
    return d
  end
end

"""
`remove_displacement(A::AbstractOperator)`

Removes the displacement of the operator.

"""
remove_displacement(A::AbstractOperator) = A

import Base: convert
function convert(::Type{T}, dom::Type, dim_in::Tuple, L::T) where {T <: AbstractOperator}
	domainType(L) != dom && error("cannot convert operator with domain $(domainType(L)) to operator with domain $dom ")
	size(L,1) != dim_in && error("cannot convert operator with size $(size(L,1)) to operator with domain $dim_in ")
	return L
end

#printing
function Base.show(io::IO, L::AbstractOperator)
	print(io, fun_name(L)*" "*fun_space(L) )
end

function fun_space(L::AbstractOperator)
	dom = fun_dom(L,2)
	codom = fun_dom(L,1)
	return dom*"->"*codom
end

function fun_dom(L::AbstractOperator,n::Int)
	dm = n == 2 ? domainType(L) : codomainType(L)
	sz = size(L,n)
	return string_dom(dm,sz)
end

function string_dom(dm::Type,sz::Tuple)
	dm_st = dm <: Complex ? " ℂ" : " ℝ"
	sz_st = length(sz) == 1 ? "$(sz[1]) " : "$sz "
	return dm_st*"^"*sz_st
end

function string_dom(dm::Tuple,sz::Tuple)
	s = string_dom.(dm,sz)
	length(s) > 3 ? s[1]*"..."*s[end]  : *(s...)
end


text:
@testset "Row" begin

@testset "Constructors and macros" begin
    @test @inferred(Row{(:A,)}((1,))) == Row{(:A,), Tuple{Int64}}((1,))

    @test @Row(A=1, B=2.0) == Row{(:A,:B), Tuple{Int,Float64}}((1,2.0))
    @test (@Row A=1 B=2.0) == Row{(:A,:B), Tuple{Int,Float64}}((1,2.0))
    @test @Row(A::Int=1, B::Float64=2.0) == Row{(:A,:B), Tuple{Int,Float64}}((1,2.0))
    @test (@Row A::Int=1 B::Float64=2.0) == Row{(:A,:B), Tuple{Int,Float64}}((1,2.0))
end

@testset "Introspection" begin
    r = Row{(:A,:B)}((1,2.0))
    @test names(r) == (:A,:B)
    @test names(Row{(:A,:B)}) == (:A,:B)
    @test names(typeof(r)) == (:A,:B)
    @test eltypes(r) == Tuple{Int64,Float64}
    @test eltypes(typeof(r)) == Tuple{Int64,Float64}
    @test length(r) == 1
    @test endof(@Row(A::Int64=1, B::Float64=2.0)) == 1
    @test ncol(r) == 2
    @test ncol(Row{(:A,:B)}) == 2
    @test ncol(typeof(r)) == 2
    @test nrow(r) == 1

    @test (show(@Row(A::Int64=1, B::Float64=2.0));println();true)

    @test names(@inferred(rename(r, Val{(:A_new, :B_new)}))) == (:A_new, :B_new)
    @test names(@inferred(rename(r, Val{:A}, Val{:A_new}))) == (:A_new, :B)

    @test @Row(A::Int64=1, B::Float64=2.0) == @Row(A::Int64=1, B::Float64=2.0)
    @test @Row(A::Int64=1, B::Float64=2.0) == @Row(B::Float64=2.0, A::Int64=1)
    @test @Row(A::Int64=1, B::Float64=2.0) != @Row(B::Float64=2.0, A::Int64=3)
end

@testset "Accessing and iterating" begin
    r = Row{(:A,:B)}((1, 2.0))
    @test copy(r) == r

    @test r[] == r
    @test r[1] == r
    @test r[:] == r

    @test r[Val{:A}] == 1
    @test r[Val{(:A,:B)}] == r
    @test r[Val{(:B,:A)}] == Row{(:B,:A)}((2.0, 1))

    @test permutecols(r, Val{(:B,:A)}) == Row{(:B,:A)}((2.0, 1))
end

@testset "Composing and Concatenating" begin
    cell_a = @Cell(A::Int64=1)
    cell_b = @Cell(B::Float64=2.0)
    cell_c = @Cell(C::Bool=true)
    row_a = @Row(A::Int64=1)
    row_b = @Row(B::Float64=2.0)

    @test @inferred(hcat(cell_a)) == Row{(:A,),Tuple{Int64}}((1,))
    @test @inferred(hcat(row_a)) == Row{(:A,),Tuple{Int64}}((1,))

    @test @inferred(hcat(cell_a, cell_b)) == Row{(:A,:B),Tuple{Int64,Float64}}((1,2.0))
    @test @inferred(hcat(cell_a, row_b)) == Row{(:A,:B),Tuple{Int64,Float64}}((1,2.0))
    @test @inferred(hcat(row_a, cell_b)) == Row{(:A,:B),Tuple{Int64,Float64}}((1,2.0))
    @test @inferred(hcat(row_a, row_b)) == Row{(:A,:B),Tuple{Int64,Float64}}((1,2.0))

    @test @inferred(hcat(cell_a, cell_b, cell_c)) == Row{(:A,:B,:C),Tuple{Int64,Float64,Bool}}((1,2.0,true))
end

end


text:
function preview_course(c::Course)
    kml = course_kml(c)
    filename = id(c) * "_preview.kml"
    open(filename,"w") do f
        for line in kml
            println(f, line)
        end
    end
end

function course_kml(c::Course)
    l = ["""<?xml version="1.0" encoding="UTF-8"?>"""]
    push!(l, """<kml xmlns="http://www.opengis.net/kml/2.2">""")
    push!(l, " <Document>")

    push!(l, """  <Style id="tee_style">""")
    push!(l, "   <IconStyle>")
    push!(l, "    <Icon>")
    push!(l, "     <href>http://maps.google.com/mapfiles/kml/paddle/grn-blank.png</href>")
    push!(l, "    </Icon>")
    push!(l, "   </IconStyle>")
    push!(l, "  </Style>")

    push!(l, """  <Style id="pin_style">""")
    push!(l, "   <IconStyle>")
    push!(l, "    <Icon>")
    push!(l, "     <href>http://maps.google.com/mapfiles/kml/paddle/blu-blank.png</href>")
    push!(l, "    </Icon>")
    push!(l, "   </IconStyle>")
    push!(l, "  </Style>")

    push!(l, """  <Style id="par_3_style">""")
    push!(l, "   <LineStyle>")
    push!(l, "    <width>1</width>")
    push!(l, "   </LineStyle>")
    push!(l, "  </Style>")

    push!(l, """  <Style id="par_not3_style">""")
    push!(l, "   <LineStyle>")
    push!(l, "    <color>ff00ffff</color>")
    push!(l, "    <width>1.5</width>")
    push!(l, "   </LineStyle>")
    push!(l, "  </Style>")

    push!(l, "  <Placemark>")
    push!(l, "   <name>$(name(c))</name>")
    push!(l, "   <Point>")
    push!(l, "    <coordinates>$(lon(c)),$(lat(c)),0</coordinates>")
    push!(l, "   </Point>")
    push!(l, "  </Placemark>")

    #Lets draw all the teepads and pins, and lines that connect them all!
    for h in holes(c)
        for p in pins(h)
            push!(l, "  <Placemark>")
            push!(l, "   <name>$(id(h) * id(p))</name>")
            push!(l, "   <styleUrl>#pin_style</styleUrl>")
            push!(l, "   <Point>")
            push!(l, "    <coordinates>$(lon(p)),$(lat(p)),0</coordinates>")
            push!(l, "   </Point>")
            push!(l, "  </Placemark>")
        end

        for t in tees(h)
            push!(l, "  <Placemark>")
            push!(l, "   <name>$(id(h) * id(t))</name>")
            push!(l, "   <styleUrl>#tee_style</styleUrl>")
            push!(l, "   <Point>")
            push!(l, "    <coordinates>$(lon(t)),$(lat(t)),0</coordinates>")
            push!(l, "   </Point>")
            push!(l, "  </Placemark>")
        end

        for  t in tees(h), p in pins(h)

            push!(l, "  <Placemark>")
            push!(l, "   <name> Hole $(id(h)), $(id(t)) ->  $(id(p))</name>")

            if par(h,t,p) == 3
                push!(l, "   <styleUrl>#par_3_style</styleUrl>")
            else
                push!(l, "   <styleUrl>#par_not3_style</styleUrl>")
            end

            push!(l, "   <LineString>")
            push!(l, "    <coordinates>")
            push!(l, "     $(lon(t)),$(lat(t))")
            push!(l, "     $(lon(p)),$(lat(p))")
            push!(l, "    </coordinates>")
            push!(l, "   </LineString>")
            push!(l, "  </Placemark>")
        end

    end

    push!(l, " </Document>")
    push!(l, "</kml>")
    l
end

function illustrate(r::Round; out="test_vis.kml")

    #write the .kml file
    xdoc = XMLDocument()

    xroot = create_root(xdoc, "kml")
    set_attribute(xroot, "xmlns", "http://www.opengis.net/kml/2.2")

    xc = new_child(xroot, "Folder")
    set_attribute(xc, "id", "Visualization")

    # Define some icon styles we will be using
    styles = [
        ("tee_style", "http://maps.google.com/mapfiles/kml/paddle/grn-blank.png", "1", "1"),
        ("pin_style", "http://maps.google.com/mapfiles/kml/paddle/blu-blank.png", "1", "1")
        ]
    for r in eachrow(DISCS_DB)
        s = (r.my_id, r.image, "1.5", "0")
        push!(styles, s)
    end
    for s in styles
        s1 = new_child(xc, "Style")
        set_attribute(s1, "id", s[1])
        s1_1 = new_child(s1, "IconStyle")
        temp_lss = new_child(s1_1, "scale")
        add_text(temp_lss, s[3])
        s1_2 = new_child(s1_1, "Icon")
        s1_3 = new_child(s1_2, "href")
        add_text(s1_3, s[2])
        temp_ls = new_child(s1, "LabelStyle")
        temp_lss = new_child(temp_ls, "scale")
        add_text(temp_lss, s[4])
    end

    # Define some line styles we will be using
    styles = [
        ("under_style", "FF00FF14", "3"),
        ("par_style", "FFFF7800", "3"),
        ("over_style", "FF1400FF", "3")
        ]
    for s in styles
        s1 = new_child(xc, "Style")
        set_attribute(s1, "id", s[1])
        s1_1 = new_child(s1, "LineStyle")
        s1_2 = new_child(s1_1, "color")
        add_text(s1_2, s[2])
        s1_3 = new_child(s1_1, "width")
        add_text(s1_3, s[3])
    end

    # create the first child, a plot of the course
    c = course(r)

    xc1 = new_child(xc, "Document")
    temp_n = new_child(xc1, "name")
    add_text(temp_n, "Course")

    # And add all the holes as children beneath that doc
    for h in holes(c)
        temp = new_child(xc1,"Document")
        temp_n = new_child(temp, "name")
        add_text(temp_n, id(h))

        for t in tees(h)
            temp_p = new_child(temp, "Placemark")
            temp_n = new_child(temp_p, "name")
            add_text(temp_n, id(h) * id(t))
            temp_s = new_child(temp_p, "styleUrl")
            add_text(temp_s, "#tee_style")
            temp_pnt = new_child(temp_p, "Point")
            temp_c = new_child(temp_pnt, "coordinates")
            add_text(temp_c, "$(lon(t)),$(lat(t)),0")
        end

        for p in pins(h)
            temp_p = new_child(temp, "Placemark")
            temp_n = new_child(temp_p, "name")
            add_text(temp_n, id(h) * id(p))
            temp_s = new_child(temp_p, "styleUrl")
            add_text(temp_s, "#pin_style")
            temp_pnt = new_child(temp_p, "Point")
            temp_c = new_child(temp_pnt, "coordinates")
            add_text(temp_c, "$(lon(p)),$(lat(p)),0")
        end

        for t in tees(h), p in pins(h)
            temp_p = new_child(temp, "Placemark")
            temp_n = new_child(temp_p, "name")
            add_text(temp_n, "$(id(h)), $(id(t)) -> $(id(p))")
            temp_ls = new_child(temp_p, "LineString")
            temp_c = new_child(temp_ls, "coordinates")
            add_text(temp_c, "$(lon(p)),$(lat(p))\n$(lon(t)),$(lat(t))")
        end

    end

    xc2 = new_child(xc, "Document")
    temp_n = new_child(xc2, "name")
    add_text(temp_n, "Round")

    for h in holes(r)
        temp = new_child(xc2,"Document")
        temp_n = new_child(temp, "name")
        add_text(temp_n, id(h.hole))

        for (i,t) in enumerate(h.throws)
            temp_p = new_child(temp, "Placemark")
            temp_n = new_child(temp_p, "name")
            add_text(temp_n, "Throw $i")
            temp_s = new_child(temp_p, "styleUrl")
            if h.res < 0
                add_text(temp_s, "#under_style")
            elseif h.res == 0
                add_text(temp_s, "#par_style")
            else
                add_text(temp_s, "#over_style")
            end
            temp_ls = new_child(temp_p, "LineString")
            temp_c = new_child(temp_ls, "coordinates")
            add_text(temp_c, "$(lon(t.loc1)),$(lat(t.loc1))\n$(lon(t.loc2)),$(lat(t.loc2))")
        end

        for (i,t) in enumerate(h.throws)
            temp_p = new_child(temp, "Placemark")
            temp_n = new_child(temp_p, "name")
            add_text(temp_n, "Throw $i")
            temp_s = new_child(temp_p, "styleUrl")
            add_text(temp_s, "#$(t.disc)")
            temp_ls = new_child(temp_p, "Point")
            temp_c = new_child(temp_ls, "coordinates")
            add_text(temp_c, "$(lon(t.loc1)),$(lat(t.loc1)),0")
        end
    end

    save_file(xdoc, out)
end


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

mutable struct ResourceGroupPatchable <: SwaggerModel
    name::Any # spec type: Union{ Nothing, String } # spec name: name
    properties::Any # spec type: Union{ Nothing, ResourceGroupProperties } # spec name: properties
    managedBy::Any # spec type: Union{ Nothing, String } # spec name: managedBy
    tags::Any # spec type: Union{ Nothing, Dict{String, String} } # spec name: tags

    function ResourceGroupPatchable(;name=nothing, properties=nothing, managedBy=nothing, tags=nothing)
        o = new()
        validate_property(ResourceGroupPatchable, Symbol("name"), name)
        setfield!(o, Symbol("name"), name)
        validate_property(ResourceGroupPatchable, Symbol("properties"), properties)
        setfield!(o, Symbol("properties"), properties)
        validate_property(ResourceGroupPatchable, Symbol("managedBy"), managedBy)
        setfield!(o, Symbol("managedBy"), managedBy)
        validate_property(ResourceGroupPatchable, Symbol("tags"), tags)
        setfield!(o, Symbol("tags"), tags)
        o
    end
end # type ResourceGroupPatchable

const _property_map_ResourceGroupPatchable = Dict{Symbol,Symbol}(Symbol("name")=>Symbol("name"), Symbol("properties")=>Symbol("properties"), Symbol("managedBy")=>Symbol("managedBy"), Symbol("tags")=>Symbol("tags"))
const _property_types_ResourceGroupPatchable = Dict{Symbol,String}(Symbol("name")=>"String", Symbol("properties")=>"ResourceGroupProperties", Symbol("managedBy")=>"String", Symbol("tags")=>"Dict{String, String}")
Base.propertynames(::Type{ ResourceGroupPatchable }) = collect(keys(_property_map_ResourceGroupPatchable))
Swagger.property_type(::Type{ ResourceGroupPatchable }, name::Symbol) = Union{Nothing,eval(Meta.parse(_property_types_ResourceGroupPatchable[name]))}
Swagger.field_name(::Type{ ResourceGroupPatchable }, property_name::Symbol) =  _property_map_ResourceGroupPatchable[property_name]

function check_required(o::ResourceGroupPatchable)
    true
end

function validate_property(::Type{ ResourceGroupPatchable }, name::Symbol, val)
end


text:
ϕ_vortex_33 =  generate_complex_potential_vortex(; pos = complex(0.0, 1.0)m, vorticity = 1.0m²/s / 2π)
ϕ_source_33 = generate_complex_potential_source(; pos = complex(3.0, 0.0)m, massflowout = 1.0m²/s)
ϕ_sink_33 = generate_complex_potential_source(; pos = -complex(3.0, 0.0)m, massflowout = -1.0m²/s)
ϕ_33(p) = ϕ_vortex_33(p) + ϕ_source_33(p) + ϕ_sink_33(p)


text:
using Test
using NQCDynamics
using NQCDynamics: Calculators
using Unitful

atoms = Atoms([:C, :H])
classical_sim = Simulation(atoms, Harmonic())
ring_polymer_sim = RingPolymerSimulation(atoms, Harmonic(), 10; temperature=100u"K")
sims = [classical_sim, ring_polymer_sim]

name(sim) = typeof(sim).name.wrapper

@testset "classical_potential_energy : $(name(sim))" for sim in sims
    u = DynamicsVariables(sim, randn(size(sim)), randn(size(sim)))
    r = DynamicsUtils.get_positions(u)
    @test DynamicsUtils.classical_potential_energy(sim, r) == DynamicsUtils.classical_potential_energy(sim, u)
end

@testset "classical_kinetic_energy : $(name(sim))" for sim in sims
    u = DynamicsVariables(sim, randn(size(sim)), randn(size(sim)))
    v = DynamicsUtils.get_velocities(u)
    @test DynamicsUtils.classical_kinetic_energy(sim, v) == DynamicsUtils.classical_kinetic_energy(sim, u)
end

@testset "divide_by_mass! : $(name(sim))" for sim in sims
    dv = randn(size(sim))
    true_value = copy(dv)
    for I in CartesianIndices(dv)
        true_value[I] /= masses(sim)[I[2]]
    end
    DynamicsUtils.divide_by_mass!(dv, masses(sim))
    @test true_value ≈ dv
end

rs = (fill(0.0, 1, 1), fill(0.0, 1, 1, 1))
calculators = (
    Calculators.Calculator(DoubleWell(), 1, Float64),
    Calculators.Calculator(DoubleWell(), 1, 1, Float64),
    )
@testset "transform_density!" for (r, calc) ∈ zip(rs, calculators)
    density = [1.0 0; 0 0]
    DynamicsUtils.transform_density!(density, calc, r, :to_adiabatic)
    @test density ≈ [0.5 0.5; 0.5 0.5]
    DynamicsUtils.transform_density!(density, calc, r, :to_diabatic)
    @test density ≈ [1.0 0; 0 0]
    @test_throws ArgumentError DynamicsUtils.transform_density!(density, calc, r, :blah)
end

@testset "initialise_adiabatic_density_matrix" for (r, calc) ∈ zip(rs, calculators)
    @testset "Adiabatic" begin
        electronics = PureState(1, Adiabatic())
        density = DynamicsUtils.initialise_adiabatic_density_matrix(electronics, calc, r)
        @test density[1,1] == 1
    end

    @testset "Diabatic" begin
        electronics = PureState(1, Diabatic())
        density = DynamicsUtils.initialise_adiabatic_density_matrix(electronics, calc, r)
        @test all(density .≈ 0.5)
    end
end


text:
# TODO What is this?
#@everywhere include(joinpath(dirname(@__FILE__),"..","src","jl"))

@everywhere using Polynomials, IterTools, Plots

srand(9)
N = 20
ran = cumsum(0.1*randn(10*N))

t = range(0, stop=1, length=N)

g6(t) = ran[floor(Int64, t*10(N-1)+1)]
ℓ = compute_transition_costs(g6, t);

Λ_0 = [create_new_pwq(EllZeroTrendFiltering.minimize_wrt_x2(ℓ[i, N], QuadraticPolynomial{Float64}(0.,0.,0.))) for i in 1:N-1];
Λ, t2, _, _, _ = @timed pwq_dp_constrained(Λ_0, ℓ, 18);


text:
@testset "455.assign-cookies.jl" begin
    @test findContentChildren([1, 2, 3], [1, 1]) == 1
    @test findContentChildren([1, 2], [1, 2, 3]) == 2
end


text:
using GridGoalMDPs
using POMDPSimulators
using POMDPPolicies
using POMDPModelTools

m = TwoAgentMDP()
p = FunctionPolicy(s->(Pos(1,0),Pos(1,0)))
for step in stepthrough(m, p, max_steps=10)
    @show step.s
    render(m, step)
end


text:
#'---
#'title : Binning Layer EAR Model Evaluation
#'author : John Waczak
#'date : 5-4-2021
#'---

#'Summon the needed Packages as well as functions from `trainModel.jl`
using Flux
using DataFrames, CSV
using Binning
using Random:shuffle!
using BSON: @save, @load
using Plots
include("./trainModel.jl")

#'Grab the data
dataPath = "../data/Oz.csv"
(Xtrain, Ytrain) = getData(dataPath)

#'Load the model
@load "./Oz-trainedWhole.bson" Model
#' Apply the model to the entire dataset
Y_model = batchModelApply(Xtrain, Model)
Y_truth = Ytrain

#'Create a scatter plot of the result
p = plot(Y_model, Y_truth, seriestype = :scatter, color=:indigo, xlabel="model", ylabel="truth", title="Avg EAR", label="")
plot!(p, [-10, 10], [-10, 10], color=:red, linestyle=:dash, label="1:1")
xlims!(p, min(Y_model...), max(Y_model...))
ylims!(p, min(Y_truth...), max(Y_truth...))

#'Save the resulting figure
savefig(p, "../figures/Avg-EAR.png")

#print out the bins
println(Model[1].b)


text:
println("Testing pkgExample ...")

using Test, pkgExample

@test pkgExample.plusone(17) == 18


text:
## options.jl : stores estimator options
#
# A type that stores estimator options shared for all Estimators.
#
# This file is part of MultilevelEstimators.jl - A Julia toolbox for
# Multilevel Monte Carlo Methods (c) Pieterjan Robbe, 2019

# valid options
default_options(::AbstractIndexSet, ::AbstractSampleMethod) = [:nb_of_warm_up_samples,
                                                               :nb_of_qoi,
                                                               :qoi_with_max_var,
                                                               :nb_of_tols,
                                                               :continuation_mul_factor,
                                                               :continuate,
                                                               :save,
                                                               :save_samples,
                                                               :verbose,
                                                               :folder,
                                                               :name,
                                                               :cost_model,
                                                               :nb_of_workers,
                                                               :nb_of_uncertainties,
                                                               :max_index_set_param,
															   :min_index_set_param]
get_valid_options(index_set, sample_method) = push!(default_options(index_set, sample_method), get_valid_options(index_set)..., get_valid_options(sample_method)...)

# additional valid options for IndexSet
default_options(::AbstractIndexSet) = [:min_splitting,
                                       :max_splitting,
                                       :robustify_bias_estimate,
                                       :do_mse_splitting,
                                       :do_regression]
get_valid_options(index_set::AbstractIndexSet) = default_options(index_set)
get_valid_options(::SL) = Symbol[]
get_valid_options(index_set::AD) = push!(default_options(index_set), :max_search_space, :penalization, :acceptance_rate)
get_valid_options(::U) = [:max_search_space, :sample_mul_factor]

# additional valid options for SampleMethod
get_valid_options(::AbstractSampleMethod) = Vector{Symbol}(undef, 0)
get_valid_options(::QMC) = [:nb_of_shifts,
                            :point_generator,
                            :sample_mul_factor]

# easy access to options fields
getindex(estimator::Estimator, s::Symbol) = estimator.options[s]


text:
# intro about metapopulations and occupancy dynamics

using MetacommunityDynamics
using DynamicGrids
using DynamicGridsGtk
using Dispersal
using Distributions
using ColorSchemes
using Plots

levinsmodel =
    RandomExtinction{:O}(0.1) +
    LevinsColonization{:O}(0.3)

latticesize = 100
initocc = 0.2
init = fill(rand() < initocc, latticesize,latticesize)

arrayout = ArrayOutput((C=initconsumer, R=initresource ), tspan=1:500)
@time sim!(arrayout, model)


text:
"System status"
function systemstatus()
  io = IOBuffer()
  versioninfo(io;verbose=true)
  vinfo = String(take!(io))
  close(io)

  status = (
    gitlog=read(`git log -1`,String),
    versioninfo=vinfo,
    datetime = now()
  )
  return status
end

project_symmetric(A) = (A+A')/2

"""
$SIGNATURES

Index in the at most `n` most frequent unique elements
of `x`. If the final number of unique elements is equal
to length(x), return `nothing`

Used in precomputation of matrix exponentials.
"""
function unique_index(x,n=length(x))
  n>0 || return nothing
  xu = unique(x)
  c = map(u->count(x.==u),xu)
  p = sortperm(c,rev=true)
  xus = xu[p]
  n = min(length(xus),n)
  xus = xus[1:n]
  length(xus)<length(x) || return nothing
  xus = xus[1:n]
  idx = indexin(x,xus)
  (idx=idx,x=xus)
end

"Is this code run from Visual Studio code?"
invscode() = occursin("code",ENV["JULIA_EDITOR"])

"""
    logpdf2pdf(lp) -> p

Convert logpdf to pdf.
"""
function logpdf2pdf(lp)
  u = exp.(lp.-maximum(lp))
  u./mean(u)
end

"
$SIGNATURES

Interpolate signal `t,x` at times `tr`.
"
function interpolate(t,y,tr;method=:nearest,extra=false)
  N = length(t)
  Nr = length(tr)
  yr = similar(y,Nr)
  issorted(t) || error("t not be sorted")
  i1 = 1; t1 = t[i1]; y1 = y[i1]
  i2 = 2; t2 = t[i2]; y2 = y[i2]
  for (i,τ) in enumerate(tr)
    if !(t[1]≤τ≤t[end])
      if !extra error("No extrapolation") end
      yr[i] = τ≤t[1] ? y[1] : y[end]
      continue
    end
    if τ>t2
      while t[i2]<τ i2+=1 end
      t2 = t[i2]; y2 = y[i2]
      i1 = i2-1
      while t[i1]>τ i1-=1 end
      t1 = t[i1]; y1 = y[i1]
    end
    yc = if method==:nearest
      t2-τ<τ-t1 ? y2 : y1
    elseif method==:linear
      (τ-t1)/(t2-t1)*(y2-y1)+y1
    else
      error("Unknown interpolation method $method")
    end
    yr[i] = yc
  end
  return yr
end

"
$SIGNATURES

Interpolate signal `t,x` at sampling interval `Tr`. Returns `tr,yr`.
"
function interpolate(t,y,Tr::Number;method=:nearest)
  tr = t[1]:Tr:t[end]
  yr = interpolate(t,y,tr;method=method)
  return tr,yr
end

"""
    randtimes([rng=GLOBAL_RNG,]Tav,n) -> t

Randomly distributed times
"""
function randtimes(rng::AbstractRNG,Tav,n,Tmin=0)
  Tavc = Tav-Tmin
  d = Exponential(Tavc)
  Δt = rand(rng,d,n).+Tmin
  t = cumsum(Δt)
end

randtimes(Tav,n) = randtimes(Random.GLOBAL_RNG,Tav,n)

"""
    semiposdef_mvrand(rng,μ,Σ)

Random draw from a multi-variate normal distribution with
a positive definite or positive semi-definite covariance matrix Σ
"""
function semiposdef_mvrand(rng,μ,Σ,tol=1e-10)
  r = if iszero(Σ)
    μ
  elseif isposdef(Σ)
    rand(rng,MvNormal(μ,Matrix(Σ)))
  else
    Qe = eigen(Σ)
    ϵ = map(Qe.values) do x
      if x≤0
          x >-tol || error("Σ has negative eigenvalue")
          0
      else
        √x*randn(rng)
      end
    end
    Qe.vectors*ϵ .+ μ
  end
  return r
end

"""
$SIGNATURES

Test process
"""
function testprocess(l)
  ## SDE roots based test processes
  r = if l==:A
    ComplexF64[0.2]
  elseif l==:B
    combineroots(Float64[],[-1/10 + 2π*0.28im])
  elseif l==:C
    combineroots([-1/18],[-1/30+2π*0.32im])
  elseif l==:D
    combineroots(Float64[],[-1/60 + 2π*0.11im,-1/22 + 2π*0.14im])
  elseif l==:E
    rr = [-1/18,-1/22,-1/140]
    ri = [-1/22 + 2π*0.14im,-1/60 + 2π*0.11im,-1/100 + 2π*0.04im]
    combineroots(rr,ri)
  elseif l==:ICMLA_A
    combineroots([-1/200],ComplexF64[])
  elseif l==:ICMLA_C
    combineroots(Float64[],[-0.10 + 2π*0.25im,-0.5 + 2π*1.5im])
  end

  if !isnothing(r)
    σy = l==:B ? 3.5 : 1
    proc = SdeProcess(SdeRoots(r;σy=σy))
    return proc
  end

  ## Other processes
  proc = if l==:ICMLA_B
    SquaredExp(0.3,1.0,0.01)
  elseif l==:ICMLA_W
    WhiteNoise()
  end

  if isnothing(proc) error("Unknown test process label $l") end

  return proc
end

"
Generate `N` observations of a discrete-time AR(1) process with parameter `a`

y[n] = a*y[n-1] + ϵ[n]

with unit standard deviation. This can be used in stress-testing SDE estimation
by generating data with negative parameter `a`, which results in an
oscillating signal that is not well described by an SDE(1) model.
"
function ar1data(rng,a,N)
  σϵ = √(1-a^2)
  y = zeros(N)
  y[1] = rand(rng,Normal())
  for n in 2:N
    μ = a*y[n-1]
    y[n] = rand(rng,Normal(μ,σϵ))
  end
  return y
end

function MCMCDiagnostics.potential_scale_reduction(θss::AbstractVector{<:NamedTuple}...)
  Tθ = eltype(first(θss))
  e1 = first(first(θss))
  Rhats = []
  for (n,T) in zip(fieldnames(Tθ),fieldtypes(Tθ))
    rh = if T <: Number
      chains = [getfield.(θs,n) for θs in θss]
      potential_scale_reduction(chains...)
    elseif T <: AbstractArray
      e1p = getfield(e1,n)
      rhc = similar(e1p)
      for i in CartesianIndices(e1p)
        chains = [getproperty.(θs,n)[i] for θs in θss]
        rhc[i] = potential_scale_reduction(chains...)
      end
      rhc
    else error("Unsupported type $T for field $n")
    end
    push!(Rhats,n=>rh)
  end
  return (;Rhats...)
end