text:
VERSION >= v"0.4.0-dev+6521" && __precompile__()

module LaTeXStrings
export LaTeXString, latexstring, @L_str, @L_mstr

using Compat
import Compat.String

# IJulia supports LaTeX output for any object with a text/latex
# writemime method, but these are annoying to type as string literals
# in Julia because of all the escaping required, e.g. "\$\\alpha +
# \\beta\$".  To simplify this, we add a new string type with a macro
# constructor, so that one can simply do L"$\alpha + \beta$".

immutable LaTeXString <: AbstractString
    s::String
end

# coercing constructor:
function latexstring(s::String)
    # the only point of using LaTeXString to represent equations, since
    # IJulia doesn't support LaTeX output other than equations, so add $'s
    # around the string if they aren't there (ignoring \$)
    return ismatch(r"[^\\]\$|^\$", s) ?
        LaTeXString(s) :  LaTeXString(string("\$", s, "\$"))
end
latexstring(s::AbstractString) = latexstring(String(s))
latexstring(args...) = latexstring(string(args...))

macro L_str(s, flags...) latexstring(s) end
macro L_mstr(s, flags...) latexstring(s) end

import Base: write, endof, getindex, sizeof, search, rsearch, isvalid, next, length, bytestring, IOBuffer, pointer
@compat import Base.show

write(io::IO, s::LaTeXString) = write(io, s.s)
@compat show(io::IO, ::MIME"application/x-latex", s::LaTeXString) = write(io, s)
@compat show(io::IO, ::MIME"text/latex", s::LaTeXString) = write(io, s)

function show(io::IO, s::LaTeXString)
    print(io, "L")
    Base.print_quoted_literal(io, s.s)
end

bytestring(s::LaTeXString) = bytestring(s.s)
endof(s::LaTeXString) = endof(s.s)
next(s::LaTeXString, i::Int) = next(s.s, i)
length(s::LaTeXString) = length(s.s)
getindex(s::LaTeXString, i::Int) = getindex(s.s, i)
getindex(s::LaTeXString, i::Integer) = getindex(s.s, i)
getindex(s::LaTeXString, i::Real) = getindex(s.s, i)
getindex(s::LaTeXString, i::UnitRange{Int}) = getindex(s.s, i)
getindex{T<:Integer}(s::LaTeXString, i::UnitRange{T}) = getindex(s.s, i)
getindex(s::LaTeXString, i::AbstractVector) = getindex(s.s, i)
sizeof(s::LaTeXString) = sizeof(s.s)
search(s::LaTeXString, c::Char, i::Integer) = search(s.s, c, i)
rsearch(s::LaTeXString, c::Char, i::Integer) = rsearch(s.s, c, i)
isvalid(s::LaTeXString, i::Integer) = isvalid(s.s, i)
pointer(s::LaTeXString) = pointer(s.s)
IOBuffer(s::LaTeXString) = IOBuffer(s.s)

# conversion to pass LaTeXString to ccall arguments
if VERSION >= v"0.4.0-dev+3710"
    import Base.unsafe_convert
else
    import Base.convert
    const unsafe_convert = Base.convert
end
@compat unsafe_convert(T::Union{Type{Ptr{UInt8}},Type{Ptr{Int8}}}, s::LaTeXString) = convert(T, s.s)
if VERSION >= v"0.4.0-dev+4603"
    unsafe_convert(::Type{Cstring}, s::LaTeXString) = unsafe_convert(Cstring, s.s)
end

end # module


text:
using Plots
using Random

Random.seed!(123)

# Solver options
tf = 5.0
N = 101
integration = :rk4
opts = SolverOptions()
opts.verbose = true
opts.cost_tolerance = 1e-8
opts.cost_tolerance_intermediate = 1e-8
opts.constraint_tolerance = 1e-8
opts.constraint_tolerance_intermediate = sqrt(opts.constraint_tolerance)
opts.gradient_norm_tolerance = 1e-12
opts.gradient_norm_tolerance_intermediate = 1e-12
opts.use_gradient_aula = false
opts.active_constraint_tolerance = 0.0
opts.penalty_scaling = 10.0
opts.penalty_initial = 1.0
opts.constraint_decrease_ratio = .25
opts.iterations = 1000
opts.iterations_outerloop = 50
opts.iterations_innerloop = 300
opts.outer_loop_update_type = :feedback

# u_min = -2
# u_max = 2
# model, obj = TrajectoryOptimization.Dynamics.pendulum!
#
# obj_con = TrajectoryOptimization.ConstrainedObjective(obj,u_min=u_min,u_max=u_max)
#
# solver_uncon = Solver(model,obj,integration=integration,N=N,opts=opts)
# solver_con = Solver(model,obj_con,integration=integration,N=N,opts=opts)

model, = TrajectoryOptimization.Dynamics.dubinscar
n, m = model.n,model.m

x0 = [0.0;0.0;0.]
xf = [0.0;1.0;0.]
tf =  3.
dt = 0.01
Qf = 100.0*Diagonal(I,n)
Q = (1e-3)*Diagonal(I,n)
R = (1e-2)*Diagonal(I,m)

obj = LQRObjective(Q, R, Qf, tf, x0, xf)

x_min = [-0.25; -0.001; -Inf]
x_max = [0.25; 1.001; Inf]

obj_con_box = TrajectoryOptimization.ConstrainedObjective(obj,x_min=x_min,x_max=x_max,use_xf_equality_constraint=true)

solver_uncon  = Solver(model, obj, integration=integration, N=N, opts=opts)
solver_con = Solver(model, obj_con_box, integration=integration, N=N, opts=opts)

U0 = rand(solver_uncon.model.m, solver_uncon.N-1)
@time results_uncon, stats_uncon = solve(solver_uncon,U0)
@time results_con, stats_con = solve(solver_con,U0)

solver_con
plot(to_array(results_con.U)[:,1:solver_con.N-1]')
plot(to_array(results_con.X)')
plot(to_array(results_con.λ[1:N-1])')

plot(log.(stats_con["max_condition_number"]))
plot(log.(stats_con["c_max"]))
plot(log.(stats_con["cost"]))
plot(log.(stats_con["gradient_norm"]))

ext: jl


text:
module test_julia_int

using Test

@test UInt8(0b11111111) == 0xff
@test UInt8(0b01111111) == 0x7f
@test UInt8(0b00111111) == 0x3f
@test UInt8(0b00011111) == 0x1f
@test UInt8(0b00001111) == 0x0f
@test UInt8(0b00000111) == 0x07
@test UInt8(0b00000011) == 0x03
@test UInt8(0b00000001) == 0x01
@test UInt8(0b00000000) == 0x00

@test trailing_ones(UInt8(0b11111111)) == 8
@test trailing_ones(UInt8(0b01111111)) == 7
@test trailing_ones(UInt8(0b00000001)) == 1

@test trailing_zeros(UInt8(0b00000000)) == 8
@test trailing_zeros(UInt8(0b10000000)) == 7
@test trailing_zeros(UInt8(0b11111110)) == 1

@test leading_ones(UInt8(0b11111111)) == 8
@test leading_ones(UInt8(0b11111110)) == 7
@test leading_ones(UInt8(0b10000000)) == 1

@test leading_zeros(UInt8(0b00000000)) == 8
@test leading_zeros(UInt8(0b00000001)) == 7
@test leading_zeros(UInt8(0b01111111)) == 1

@test widen(Int64) === Int128
@test 1 isa Int64
@test typeof(widen(1)) === Int128
@test signed(UInt64) === Int64
@test unsigned(Int64) === UInt64
@test signbit(1) === false
@test signbit(-1)

bits(::Int64) = 64
bits(::Int32) = 32
@test bits(1) == 64
@test bits(Int32(1)) == 32

end # module test_julia_int

module test_julia_integer

using Test

f(n::Integer) = (:integer, n)
@test f(1) == (:integer, 1)

f(n::Int) = (:int, n)
@test f(1) == (:int, 1)
@test f(true) == (:integer, true)

end # module test_julia_integer

module test_julia_bswap

using Test

@test bswap(0x3039) == 0x3930
@test (bswap ∘ bswap)(0x3039) == 0x3039

end # module test_julia_bswap


text:
# Use baremodule to shave off a few KB from the serialized `.ji` file
baremodule Infernal_jll
using Base
using Base: UUID
import JLLWrappers

JLLWrappers.@generate_main_file_header("Infernal")
JLLWrappers.@generate_main_file("Infernal", UUID("5109d7c4-8f61-540a-bdd6-fb9d06080a8f"))
end  # module Infernal_jll


text:
using IonCLI

IonCLI.command_main(["-h"])


text:
using Test

# Identify files in test/ that are testing matching files in src/
#  - src/Source.jl will be matched by test/test_Source.jl

filebase = map(file -> replace(file, r"(.*).jl$" => s"\1"),
                filter(file -> occursin(r".*\.jl$", file),
                       readdir("../src")))
testbase = map(file -> replace(file, r"test_(.*).jl$" => s"\1"),
                filter(str -> occursin(r"^test_.*\.jl$", str), readdir()))

println()
@info "Running tests for files:"
for t in testbase
    println("    = $t.jl")
end
println()

@testset "* Testing $t.jl" for t in testbase
    include("test_$t.jl")
end

# Identify tests with no matching file
superfluous = filter(f -> f ∉ filebase, testbase)
if length(superfluous) > 0
    println()
    @info "Potentially superfluous tests:"
    for f in superfluous
        println("    + $f.jl")
    end
    println()
end

# Identify files with no matching test
missing = filter(f -> f ∉ testbase, filebase)
if length(missing) > 0
    println()
    @info "Potentially missing tests:"
    for f in missing
        println("    - $f.jl")
    end
    println()
end

# Identify files that are cross-validating results against other packages
# test/pkg_Package.jl should validate results against the Package package

pkgbase = map(file -> replace(file, r"pkg_(.*).jl$" => s"\1"),
                   filter(str -> occursin(r"^pkg_.*\.jl$", str),
                          readdir()))

if length(pkgbase) > 0
    @info "Running cross-validation against:"
    for p in pkgbase
        println("    = $p")
    end
    println()

    @testset " * Validating against $p" for p in pkgbase
        include("pkg_$p.jl")
    end
end


text:
using BenchmarkTools

function my_append()

    v = []

    for i in 1:10000
        push!(v, 5)
    end
end

function my_prealloc()

    v = Vector{Int64}(undef, 10000)

    for i in 1:10000
        v[i] = 5
    end
end

# @benchmark my_append()
@benchmark my_prealloc()

max_stars_repo_name: nunesmelo/djs-office-hours

max_issues_repo_name: nunesmelo/djs-office-hours

max_forks_repo_name: nunesmelo/djs-office-hours


text:
function part1(input)
    timestamp, bus_list = readlines(input)
    time = parse(Int, timestamp)
    buses = parse.(Int, filter(!=("x"), split(bus_list, ",")))
    wait_times = mod.(-time, buses)
    _, i = findmin(wait_times)
    return wait_times[i] * buses[i]
end

function part2(input)
    _, bus_list = readlines(input)
    buses = tryparse.(Int, split(bus_list, ","))
    modulo = BigInt(1)
    remainder = BigInt(0)
    for (i, n) in enumerate(buses)
        isnothing(n) && continue
        modulo, remainder = crt(modulo, remainder, BigInt(n), BigInt(1 - i))
    end
    return remainder
end

function crt(n1, r1, n2, r2)
    d, u, v = gcdx(n1, n2)
    n = n1 * n2 ÷ d
    return n, mod((r1 * n2 * v + r2 * n1 * u) ÷ d, n)
end


text:
#    Copyright (c) 2021 J.A. Duffek
#    Copyright (c) 2000 D.M. Spink
#
#    This program is free software; you can redistribute it and/or modify
#    it under the terms of the GNU General Public License as published by
#    the Free Software Foundation; either version 2 of the License, or
#    (at your option) any later version.
#
#    This program is distributed in the hope that it will be useful,
#    but WITHOUT ANY WARRANTY; without even the implied warranty of
#    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
#    GNU General Public License for more details.
#
#    You should have received a copy of the GNU General Public License
#    along with this program; if not, see <http://www.gnu.org/licenses/>.

using Plots;

"""
    demo4surf()

Demonstration of a bilinear surface.

# Examples:
```julia
julia> demo4surf()
```
"""
function demo4surf()
# make and draw nurbs test surface
srf = nrb4surf([0.0;0.0;0.5],[1.0;0.0;-0.5],[0.0;1.0;-0.5],[1.0;1.0;0.5]);
nrbplot(srf,[10;10],
        title = "Construction of a bilinear surface.",
        c = :winter,
        legend = nothing,
        linewidth = 0.5,
        linecolor = :black,
        camera=[-30,30]);
end # demo4surf


text:
using BinaryBuilder, Pkg

const YGGDRASIL_DIR = "../../.."
include(joinpath(YGGDRASIL_DIR, "fancy_toys.jl"))
include(joinpath(YGGDRASIL_DIR, "platforms", "cuda.jl"))

name = "CUDA_Runtime"
version = v"0.1"

cuda_versions = [v"10.0", v"10.2",
                 v"11.0", v"11.1", v"11.2", v"11.3", v"11.4", v"11.5", v"11.6"]

augment_platform_block = """
    using Base.BinaryPlatforms

    $(CUDA.augment)

    augment_cuda_toolkit!(platform) = augment_cuda_toolkit!(platform, $cuda_versions)

    function augment_platform!(platform::Platform)
        augment_cuda_toolkit!(platform)
    end"""

# determine exactly which tarballs we should build
builds = []
for cuda_version in cuda_versions
    cuda_tag = "$(cuda_version.major).$(cuda_version.minor)"
    include("build_$(cuda_tag).jl")

    for platform in platforms
        augmented_platform = deepcopy(platform)
        augmented_platform[CUDA.platform_name] = CUDA.platform(cuda_version)

        should_build_platform(triplet(augmented_platform)) || continue
        push!(builds, (;
            dependencies, script, products,
            platforms=[augmented_platform],
        ))
    end
end

# don't allow `build_tarballs` to override platform selection based on ARGS.
# we handle that ourselves by calling `should_build_platform`
non_platform_ARGS = filter(arg -> startswith(arg, "--"), ARGS)

# `--register` should only be passed to the latest `build_tarballs` invocation
non_reg_ARGS = filter(arg -> arg != "--register", non_platform_ARGS)

for (i,build) in enumerate(builds)
    build_tarballs(i == lastindex(builds) ? non_platform_ARGS : non_reg_ARGS,
                   name, version, [], build.script,
                   build.platforms, build.products, build.dependencies;
                   julia_compat="1.6", lazy_artifacts=true,
                   augment_platform_block)
end

max_stars_repo_name: MichelJuillard/Yggdrasil

max_issues_repo_name: MichelJuillard/Yggdrasil

max_forks_repo_name: MichelJuillard/Yggdrasil


text:
__precompile__(true)

module Learn

using Reexport
@reexport using LearnBase
@reexport using LossFunctions
@reexport using Transformations
@reexport using ObjectiveFunctions
@reexport using PenaltyFunctions
@reexport using StochasticOptimization
# @reexport using MLDataUtils
@reexport using MLMetrics

end # module


text:
"""
    simple_returns(prices::Vector)

Calculates the simple returns series based on the provided time series of `N` prices.
Please note that a) the provided prices series should incorporate
all corporate actions, dividends etc., and b) should be regularly spaced, for example hourly or daily data.

# Arguments
- `prices`: Vector of adjusted prices accounting for all corporate actions etc.

# Returns
Vector of simple returns with size `N-1`.

# Examples
```jldoctest
julia> prices = [100.0, 101.0, 100.5, 99.8, 101.3]
5-element Array{Float64,1}:
 100.0
 101.0
 100.5
  99.8
 101.3

julia> simple_returns(prices)
4-element Array{Float64,1}:
  0.010000000000000009
 -0.004950495049504955
 -0.006965174129353269
  0.01503006012024044
```
"""
@inline function simple_returns(prices::Vector)
    a = @view prices[2:end]
    b = @view prices[1:end-1]
    (a ./ b) .- 1.0
end

"""
    simple_returns(dates::Vector{DateTime}, prices::Vector{T}; drop_overnight::Bool=false)

Calculates the simple returns series based on the provided time series of prices.
Please note that a) the provided prices series should incorporate
all corporate actions, dividends etc., and b) should be regularly spaced, for example hourly or daily data.
If the parameter `drop_overnight` is set to `true`, overnight returns will be dropped.

# Arguments
- `dates`:              Vector of `DateTime` objects for prices.
- `prices`:             Vector of adjusted prices accounting for all corporate actions etc.
- `drop_overnight`:     Boolean value indicating whether to drop overnight returns or not (default=false).

# Returns
Tuple of type `Tuple{Vector{DateTime}, Vector{T}, Vector{Int64}}` with dates, returns and indices of original time series.

# Examples
```
julia> dates = [DateTime(2020, 1, 1, 1), DateTime(2020, 1, 1, 2), DateTime(2020, 1, 2, 1), DateTime(2020, 1, 2, 2)]
4-element Array{DateTime,1}:
 2020-01-01T01:00:00
 2020-01-01T02:00:00
 2020-01-02T01:00:00
 2020-01-02T02:00:00

julia> prices = [100.0, 101.0, 100.5, 100.75]
4-element Array{Float64,1}:
 100.0
 101.0
 100.5
 100.75

julia> simple_returns(dates, prices; drop_overnight=true)
([DateTime("2020-01-01T02:00:00"), DateTime("2020-01-02T02:00:00")], [0.010000000000000009, 0.0024875621890547706], [2, 4])
```
"""
function simple_returns(
    dates           ::Vector{DateTime},
    prices          ::Vector{T};
    drop_overnight  ::Bool=false
)::Tuple{Vector{DateTime}, Vector{T}, Vector{Int64}} where T
    N = length(prices)
    if drop_overnight
        a = @view dates[2:N]
        b = @view dates[1:N-1]
        idxs = findall(vcat(false, Dates.days.(a) .== Dates.days.(b)))
        a2 = @view prices[idxs]
        b2 = @view prices[idxs .- 1]
        return dates[idxs], (a2 ./ b2) .- 1.0, idxs
    else
        return dates[2:N], simple_returns(prices), collect(2:N)
    end
end

"""
    simple_returns(prices::Matrix)

Calculates the simple returns series for each column of the provided prices matrix, where each column denotes a prices series of an asset, e.g. a stock.

Please note that a) the provided prices series should incorporate all corporate actions,
dividends etc. and b) should be regularly spaced, for example hourly or daily data.

# Arguments
- `prices`: Matrix of size `[N x k]` of adjusted prices accounting for all corporate actions etc., where each column denotes a price series of an asset, e.g. a stock.

# Returns
Matrix of simple returns with size `[(N-1) x k]`.

# Examples
```jldoctest
julia> prices = [100.0 99.0; 101.0 99.5; 100.5 99.8; 99.8 99.2; 101.3 100.0]
5×2 Array{Float64,2}:
 100.0   99.0
 101.0   99.5
 100.5   99.8
  99.8   99.2
 101.3  100.0

julia> simple_returns(prices)
4×2 Array{Float64,2}:
  0.01         0.00505051
 -0.0049505    0.00301508
 -0.00696517  -0.00601202
  0.0150301    0.00806452
```
"""
@inline function simple_returns(prices::Matrix)
    a = @view prices[2:end, :]
    b = @view prices[1:end-1, :]
    (a ./ b) .- 1.0
end

"""
    log_returns(prices::Vector)

Calculates the log-return series based on the provided time series of `N` prices.
Please note that a) the provided prices series should incorporate
all corporate actions, dividends etc., and b) should be regularly spaced, for example hourly or daily data.

# Arguments
- `prices`: Vector of adjusted prices accounting for all corporate actions etc.

# Returns
Vector of log-returns with size `N-1`.

# Examples
```jldoctest
julia> prices = [100.0, 101.0, 100.5, 99.8, 101.3]
5-element Array{Float64,1}:
 100.0
 101.0
 100.5
  99.8
 101.3

julia> log_returns(prices)
4-element Array{Float64,1}:
  0.009950330853168092
 -0.004962789342129014
 -0.006989544181712186
  0.014918227937219366
```
"""
@inline function log_returns(prices::Vector)
    a = @view prices[2:end]
    b = @view prices[1:end-1]
    log.(a ./ b)
end

"""
    log_returns(prices::Matrix)

Calculates log-return series for each column of the provided prices matrix, where each column denotes a prices series of an asset, e.g. a stock.

Please note that a) the provided prices series should incorporate all corporate actions,
dividends etc. and b) should be regularly spaced, for example hourly or daily data.

# Arguments
- `prices`: Matrix of size `[N x k]` of adjusted prices accounting for all corporate actions etc., where each column denotes a price series of an asset, e.g. a stock.

# Returns
Matrix of log-returns with size `[(N-1) x k]`.

# Examples
```jldoctest
julia> prices = [100.0 99.0; 101.0 99.5; 100.5 99.8; 99.8 99.2; 101.3 100.0]
5×2 Array{Float64,2}:
 100.0   99.0
 101.0   99.5
 100.5   99.8
  99.8   99.2
 101.3  100.0

julia> log_returns(prices)
4×2 Array{Float64,2}:
  0.00995033   0.00503779
 -0.00496279   0.00301054
 -0.00698954  -0.00603017
  0.0149182    0.00803217
```
"""
@inline function log_returns(prices::Matrix)
    a = @view prices[2:end, :]
    b = @view prices[1:end-1, :]
    log.(a ./ b)
end


text:
# Fractional Order Chaotic examples come from
#=
@inproceedings{Petr2011FractionalOrderNS,
title={Fractional-Order Nonlinear Systems: Modeling, Analysis and Simulation},
author={Ivo Petr{\'a}{\vs}},
year={2011}
}
=#
#

#=======Fractional Order Liu System=======#
using FractionalDiffEq

h=0.005
alpha = [0.95, 0.95, 0.95]
x0 = [0.2, 0, 0.5]
tf=100
function Liu!(du, u, p, t)
    a, b, c, e, n, m = 1, 2.5, 5, 1, 4, 4
    du[1] = -a*u[1]-e*u[2]^2
    du[2] = b*u[2]-n*u[1]*u[3]
    du[3] = -c*u[3]+m*u[1]*u[2]
end
prob = FODESystem(Liu!, alpha, x0, tf)
result = solve(prob, h, GLWithMemory())

using Plots
plot3d(result[:, 1], result[:, 2], result[:, 3], title="Fractional Order Liu System")

#=======Fractional Order Duffing System=======#
using FractionalDiffEq

h=0.005
alpha = [0.9, 1]
x0 = [0.21, 0.31]
tf=100
function Duffing!(du, u, p, t)
    α, δ, ω = 0.15, 0.3, 1
    du[1] = u[2]
    du[2] = u[1]-α*u[2]-u[1]^3+δ*cos(ω*t)
end
prob = FODESystem(Duffing!, alpha, x0, tf)
result = solve(prob, h, GLWithMemory())

using Plots
plot(result[:, 1], result[:, 2], title="Fractional Order Duffing System")

#=======Fractional Order Chen System=======#
using FractionalDiffEq

h=0.005
alpha = [0.9, 0.9, 0.9]
x0 = [-9, -5, 14]
tf=100
function Chen!(du, u, p, t)
    a, b, c, d = 35, 3, 28, -7
    du[1] = a*(u[2]-u[1])
    du[2] = d*u[1]-u[1]*u[3]+c*u[2]
    du[3] = u[1]*u[2]-b*u[3]
end
prob = FODESystem(Chen!, alpha, x0, tf)
result = solve(prob, h, GLWithMemory())

using Plots
plot3d(result[:, 1], result[:, 2], result[:, 3], title="Fractional Order Chen System")

#=======Fractional Order Van der Pol Oscillator=======#
using FractionalDiffEq

h=0.005
alpha = [1.2, 0.8]
x0 = [0.2, -0.2]
tf=60
function VanderPol!(du, u, p, t)
    ϵ = 1
    du[1] = u[2]
    du[2] = -u[1]-ϵ*(u[1]^2-1)*u[2]
end
prob = FODESystem(VanderPol!, alpha, x0, tf)
result = solve(prob, h, GLWithMemory())

using Plots
plot(result[:, 1], result[:, 2], title="Fractional Order Van der Pol Oscillator")

#=======Fractional Order Lotka Volterra System=======#
using FractionalDiffEq

h=0.005
alpha = [0.95, 0.95, 0.95]
x0 = [1, 1.4, 1]
tf=50
function LotkaVolterra!(du, u, p, t)
    a, b, c, d, e, p, s = 1, 1, 1, 1, 2, 3, 2.7
    du[1] = a*u[1]+e*u[1]^2-b*u[1]*u[2]-s*u[3]*u[1]^2
    du[2] = -c*u[2]+d*u[1]*u[2]
    du[3] = -p*u[3]+s*u[3]*u[1]^2
end
prob = FODESystem(LotkaVolterra!, alpha, x0, tf)
result = solve(prob, h, GLWithMemory())

using Plots
plot3d(result[:, 1], result[:, 2], result[:, 3], title="Fractional Order Lotka Volterra System")

#=======Fractional Order Lu System=======#
using FractionalDiffEq

h=0.005
alpha = [0.985, 0.99, 0.98]
x0 = [0.2, 0.5, 0.3]
tf=60
function Lu!(du, u, p, t)
    a, b, c = 36, 3, 20
    du[1] = a*(u[2]-u[1])
    du[2] = -u[1]*u[3]+c*u[2]
    du[3] = u[1]*u[2]-b*u[3]
end
prob = FODESystem(Lu!, alpha, x0, tf)
result = solve(prob, h, GLWithMemory())

using Plots
plot3d(result[:, 1], result[:, 2], result[:, 3], title="Fractional Order Lu System")

#=======Fractional Order Rossler System=======#
using FractionalDiffEq

h=0.005
alpha = [0.9, 0.85, 0.95]
x0 = [0.5, 1.5, 0.1]
tf=120
function Rossler!(du, u, p, t)
    a, b, c = 0.5, 0.2, 10
    du[1] = -u[2]-u[3]
    du[2] = u[1]+a*u[2]
    du[3] = b+u[3]*(u[1]-c)
end
prob = FODESystem(Rossler!, alpha, x0, tf)
result = solve(prob, h, GLWithMemory())

using Plots
plot3d(result[:, 1], result[:, 2], result[:, 3], title="Fractional Order Rossler System")

#=======Fractional Order Arneodo System=======#
using FractionalDiffEq

h=0.005
alpha = [0.97, 0.97, 0.96]
x0 = [-0.2, 0.5, 0.2]
tf=100
function Arneodo!(du, u, p, t)
    b1, b2, b3, b4 = -5.5, 3.5, 0.8, -1.0
    du[1] = u[2]
    du[2] = u[3]
    du[3] = -b1*u[1]-b2*u[2]-b3*u[3]+b4*u[1]^3
end
prob = FODESystem(Arneodo!, alpha, x0, tf)
result = solve(prob, h, GLWithMemory())

using Plots
plot3d(result[:, 1], result[:, 2], result[:, 3], title="Fractional Order Arneodo System")

#=======Fractional Order Genesio-Tesi System=======#
using FractionalDiffEq

h=0.005
alpha = [0.95, 0.95, 0.95]
x0 = [-0.1, 0.5, 0.2]
tf=100
function GenesioTesi!(du, u, p, t)
    b1, b2, b3, b4 = 1, 1.1, 0.4, 1.0
    du[1] = u[2]
    du[2] = u[3]
    du[3] -b1*u[1]-b2*u[2]-b3*u[3]+b4*u[1]^2
end
prob = FODESystem(GenesioTesi!, alpha, x0, tf)
result = solve(prob, h, GLWithMemory())

using Plots
plot3d(result[:, 1], result[:, 2], result[:, 3], title="Fractional Order Genesio-Tesi System")


text:
module  CumEx

#this is useful if we have mdoules
if pwd() ∉ LOAD_PATH
  push!(LOAD_PATH,pwd())
end

using DataFrames, Distributions, StatsBase, GZip, JLD, Gadfly

#=#NOTE##############
Some companies have multiple PERMNOs. The returns are different, so I am leaving them as is.
=####################
const DATA_PATH = pwd() * "\\data" #path where we store the data
const DATA_FILE_NAME = "ReturnsDivOnly1926-2017" #data file name
const DATE_FORMAT_STR = "yyyymmdd"
const WINSOR_MAX = 0.2
const NET_OF_MARKET = true
const MIN_DIVIDEND = 10.0^-4.0 #1/100 of a penny
const MIN_PRICE = 1.0
const EXC_SIC_CODES = nothing

#useful type
const CTSym = Union{Symbol,Void}

#This function converts a string datavector to a real number
#IN: A datavector string
#OUT: A data vector with all string numbers converted to numbers,
# and everything else converted to an NA
function parseDataVector(strings::DataVector{String})::DataVector{Float64}
  N::Int = length(strings)

  vals::DataVector{Float64} = DataVector{Float64}(Vector{Float64}(N),trues(N))

  #threading makes this much faster
  @fastmath @inbounds Threads.@threads for i::Int ∈ 1:N
    parsed::Union{Symbol, Float64} = parse(strings[i])
    vals[i] = typeof(parsed)<:Float64?parsed:NA
  end

  return vals
end

#in this case, do nothing
parseDataVector(vals::DataVector{Float64})::DataVector{Float64}  = vals

#This function winsorizes a data column
#IN: A datavector string
#OUT: A data vector with all string numbers converted to numbers,
# and everything else converted to an NA
winsorize(col::DataVector{Float64}, minVal::Float64, maxVal::Float64)::DataVector{Float64} =
  min.(max.(col, minVal),maxVal)

#same as above but works at the element level
winsorize(val::T, minVal::Float64, maxVal::Float64) where T<:Union{NAtype,Float64} =
  min.(max.(val, minVal),maxVal)

#this function pre-processes the data files
#IN: A switch for reading or re-reading the csv
#OUT: nothing. A clean and compressed julia serial file has all of the data
function processStockData(;newBinary::Bool = true)::Void
  dtFormat::DateFormat = DateFormat(DATE_FORMAT_STR) ## a date format object

  if !isfile("$DATA_PATH\\$DATA_FILE_NAME.jls.gz") || newBinary

  #extract from the zip file, package into dataframe and serialize into zip
  #Saves a lot of time not to have to do this each time
    igStream::GZipStream = gzopen("$DATA_PATH\\$DATA_FILE_NAME.csv.gz")
    ogStream::GZipStream = gzopen("$DATA_PATH\\$DATA_FILE_NAME.jls.gz","w")
    serialize(ogStream, readtable(igStream))

    close(igStream)
    close(ogStream)
  end

  igStream = gzopen("$DATA_PATH\\$DATA_FILE_NAME.jls.gz")
  stockDF::DataFrame  = deserialize(igStream)
  close(igStream)

  #unfortunately, some of these values have error codes
  stockDF = stockDF[completecases(stockDF[:,[:RET,:RETX]]),:]
  stockDF[:,:RET] = parseDataVector(stockDF[:,:RET])
  stockDF[:,:RETX] = parseDataVector(stockDF[:,:RETX])
  stockDF = stockDF[completecases(stockDF[:,[:RET,:RETX,:PRC]]),:] #can get more NAs from before

  stockDF = stockDF[stockDF[:,:DIVAMT] .≥ MIN_DIVIDEND,:] #Filter out extremely small dividends
  stockDF = stockDF[stockDF[:,:PRC] .≥ MIN_PRICE,:] #Filter out penny stocks

  if EXC_SIC_CODES ≠ nothing
    stockDF = stockDF[(stockDF[:,:SICCD] .> maximum(EXC_SIC_CODES)) | (
      stockDF[:,:SICCD] .< minimum(EXC_SIC_CODES)),:] #Filter out financials
  end

  stockDF = deepcopy(stockDF[:,[:PERMNO,:date,:PERMCO, :DIVAMT,:PRC,:SHROUT, :RET,:RETX,:vwretd]])
  stockDF[:,:date] = ((x::Int)->Date("$x",dtFormat)).(stockDF[:,:date])
  N::Int = size(stockDF,1)

  #generate the derived columns
  stockDF[:,:divYield] = stockDF[:,:RET] .- stockDF[:,:RETX]
  stockDF[:,:RETXW] = winsorize(stockDF[:,:RETX],-1.0*WINSOR_MAX,WINSOR_MAX)
  stockDF[:,:divYieldW] = winsorize(stockDF[:,:divYield],-1.0*WINSOR_MAX,WINSOR_MAX)

  stockDF[:,:marketCap] = stockDF[:,:SHROUT] .* stockDF[:,:PRC]

  stockDF[:,:year] = (Dates.year).(stockDF[:,:date])
  stockDF[:,:month] = (Dates.month).(stockDF[:,:date])
  stockDF[:,:quarter] = (Dates.quarterofyear).(stockDF[:,:date])
  stockDF[:,:quarterOfYear] = (Dates.firstdayofquarter).(stockDF[:,:date])
  stockDF[:,:monthOfYear] = (Dates.firstdayofmonth).(stockDF[:,:date])

  ogStream = gzopen("$DATA_PATH\\$(DATA_FILE_NAME)_clean.jls.gz","w")
  serialize(ogStream, stockDF)
  close(ogStream)

  #showcols(stockDF)
  return nothing
end

#this function is meant to be used internally in aggDivCumEx to build the time series
#IN: a dataframe or subdataframe of input data and a weighting vector
#OUT: the weighted sum of the input vector. This is a weighted mean if Σw=1
function getNegROverDNet(df::T, weightVec::Vector{Float64})::Float64 where {T<:AbstractDataFrame}

  N::Int = size(df, 1)

  ans::Float64 = sum(-(df[:,:RETXW].-df[:,:vwretd])./df[:,:divYieldW] .* weightVec)
  ans /= sum(weightVec)
  return ans
end

function getDivLevels(df::T, weightVec::Vector{Float64})::Float64 where {T<:AbstractDataFrame}

  ans::Float64 = sum(df[:,:DIVAMT] .* weightVec)
  #ans /= sum(weightVec)

  return ans
end

#aggregates the Elton-Gruber 1970 statistic based on the supplied groupBy column and weight
#IN: The dataframe, the column to split by, weights, and an aggregation function
#OUT: A dataframe aggregated according to the input parameters
#NOTE: Default options supplied for testing purposes
function aggDivCumEx(stockDF::DataFrame, weightCol::Symbol; byCol::Symbol=:year,
  statFunc::Function = getNegROverDNet)::DataFrame

  #apply the supplied function to each subsection
  outDF::DataFrame = by(stockDF, byCol) do subDF::SubDataFrame
    DataFrame(EGStat = statFunc(subDF,Vector{Float64}(subDF[:,weightCol])))
  end

  return outDF
end

#this version allows for a vector of weighting columns.
#The results are merged into a single dataframe
#IN: the source data column, the weighting columns, the by column for aggregation, and
# the aggregation function
#OUT: The merged data frame.
#NOTE: Each weighting column can be void, and if so this creates a default equal-weighted column
function aggDivCumEx(stockDF::DataFrame, weightCols::Vector{T}=[nothing]; byCol::Symbol=:year,
  statFunc::Function = getNegROverDNet)::DataFrame where T<:CTSym

  #set unit weights if no other weights
  if nothing ∈ weightCols
    weightCol = :equalWeights
    stockDF[:,:equalWeights] = ones(size(stockDF,1))
  end

  #replace the nothing column
  weightCols =
    [weightCols[i]==nothing?(:equalWeights):(weightCols[i]) for i::Int ∈ 1:length(weightCols)]

  #Create the dataframe from the first aggregation
  outDF = aggDivCumEx(stockDF, :equalWeights, byCol=byCol, statFunc = getDivLevels)
  rename!(outDF, :EGStat, :Levels) #need to set the name for each column

  #if more than one weight is provided, execute again and append to the dataframe

  for weightCol::Symbol ∈ weightCols[1:end]
    outDF =
      [outDF aggDivCumEx(stockDF, weightCol, byCol=byCol, statFunc = getNegROverDNet)[:,:EGStat]]
    #println(outDF)
    rename!(outDF, :x1, Symbol(:EGStat,"_", weightCol))
  end

  return outDF
end

#This is a script for running the aggregation. Uses the output of processStockData
#IN: Nothing, but reads in a dataframe as a serial file
#OUT: Nothing, but writes the aggregated tables to the disk
function getCumExTablesScript()::Void
  igStream = gzopen("$DATA_PATH\\$(DATA_FILE_NAME)_clean.jls.gz")
  stockDF::DataFrame  = deserialize(igStream)
  close(igStream)

  #set the dimensions of the cum-ex aggregation
  weightCols::Vector{CTSym} = [nothing, :divYieldW, :marketCap]
  byCols::Vector{Symbol} = [:year, :quarterOfYear, :monthOfYear]
  aggTables::Vector{DataFrame} = Vector{DataFrame}()

  #get the dataframes
  for byCol::Symbol ∈ byCols
    push!(aggTables, aggDivCumEx(stockDF, weightCols, byCol=byCol, statFunc = getNegROverDNet))
  end

  #write out the dataframes
  for i::Int ∈ 1:length(aggTables)
    pth::String = "$DATA_PATH\\$(DATA_FILE_NAME)_$(byCols[i]).csv"
    println("Writing: $pth")
    writetable(pth, aggTables[i])
  end

  println("All tables writtens successfully.")

  return nothing
end

#run the appropriate functions
@time begin
processStockData(newBinary = false)
getCumExTablesScript()
end

end

max_stars_repo_name: clintonTE/CCA

max_issues_repo_name: clintonTE/CCA

max_forks_repo_name: clintonTE/CCA


text:
    # data =

    # m = 2n+1
    # E_check = TD.empirical_coefs(c_data,(m,))

    # # reconstruction:
    # τ = TD.thorin_moments(vec(c_data),t_star,m)
    # κ = τ .* factorial.([big(0),big.(0:m-1)...])
    # μ = TD.univ_μ_from_κ(κ)
    # A = TD.get_precomp(BigFloat,sum(m)).LAGUERRE .* sqrt(2)
    # E = A'μ[1:end-1]

    # # E should be equal to E_check.

    # and
    #matrix_A(m)' == A


text:
function insertionsort!(data::Array{Int64})
	for i in 2 : length(data)
		current_value = data[i]
		for j in i - 1 : -1 : 1
			if current_value < data[j]
				data[j + 1] = data[j]
				if j == 1
					data[j] = current_value
				end
			else
				data[j + 1] = current_value
				break
			end
		end
	end
end

function builtinsort!(data::Array{Int64})
	sort!(data)
end

function quicksort!(data::Array, start=-1, stop=-1, pivot=-1)
	if start == -1
		start = 1
		stop = length(data)
		pivot = 0
		for i in data
			pivot += i
		end
		pivot = pivot / stop
	end

	print(start)
	print(", ")
	print(stop)
	print(", ")
	println(pivot)

	if stop - start <= 2
		if stop - start == 2
			if data[start] > data[start + 1]
				buffer = data[start]
				data[start] = data[start + 1]
				data[start + 1] = buffer
			end
			if data[start + 1] > data[start + 2]
				buffer = data[start + 1]
				data[start + 1] = data[start + 2]
				data[start + 2] = buffer
			end
			if data[start] > data[start + 1]
				buffer = data[start]
				data[start] = data[start + 1]
				data[start + 1] = buffer
			end
			println(data)
			return
		end

		if stop - start == 1
			if data[start] > data[start + 1]
				buffer = data[start]
				data[start] = data[start + 1]
				data[start + 1] = buffer
			end
			println(data)
			return
		end

		if stop - start == 0
			return
		end
	end

	leftPointer = start
	rightPointer = stop

	done = false

	while !done
		print(leftPointer)
		print(", ")
		println(rightPointer)

		left = data[leftPointer] <= pivot
		right = data[rightPointer] > pivot

		if rightPointer - leftPointer <= 0
			done = true
		elseif !left && !right
			buffer = data[leftPointer]
			data[leftPointer] = data[rightPointer]
			data[rightPointer] = buffer
			if rightPointer - leftPointer > 1
				leftPointer += 1
				rightPointer -= 1
			else
				done = true
			end
		elseif !left && right
			rightPointer -= 1
		elseif left && !right
			leftPointer += 1
		elseif rightPointer - leftPointer > 1
			leftPointer += 1
			rightPointer -= 1
		else
			done = true
		end
		println(data)
	end

	pivotPointer = leftPointer

	if leftPointer == rightPointer && data[leftPointer] > pivot
		pivotPointer -= 1
	end

	quicksort!(data, start, pivotPointer, pivot / 2)
	quicksort!(data, pivotPointer + 1, stop, pivot * 3 / 2)
end


text:
import Pkg

function _generate_packagecompiler_run_content(config::Config)
    julia_cpu_target = config.julia_cpu_target
    make_sysimage = config.make_sysimage
    exclude_packages_from_sysimage = config.exclude_packages_from_sysimage
    lines = String[
        "import Pkg",

        "all_stdlib_uuids = collect(keys(Pkg.Types.stdlibs()))",
        "if Base.VERSION >= v\"1.8-\"",
        "all_stdlib_names = first.(collect(values(Pkg.Types.stdlibs())))",
        "else",
        "all_stdlib_names = collect(values(Pkg.Types.stdlibs()))",
        "end",

        "pkgnames = Vector{String}(undef, 0)",
        "for (uuid, info) in Pkg.dependencies()",
        "    if !(info.name in $(exclude_packages_from_sysimage))",
        "        push!(pkgnames, info.name)",
        "    end",
        "end",
        "sort!(pkgnames)",

        "pkgnames_nonstdlib = Vector{String}(undef, 0)",
        "pkgnames_stdlib = Vector{String}(undef, 0)",
        "for name in pkgnames",
        "    if name in all_stdlib_names",
        "        push!(pkgnames_stdlib, name)",
        "    else",
        "        push!(pkgnames_nonstdlib, name)",
        "    end",
        "end",

        "pkgnames_symbols = Symbol.(pkgnames)",
        "pkgnames_nonstdlib_symbols = Symbol.(pkgnames_nonstdlib)",
        "pkgnames_stdlib_symbols = Symbol.(pkgnames_stdlib)",
        "println(\"The length of pkgnames_symbols is \$(length(pkgnames_symbols))\")",
        "println(\"The length of pkgnames_nonstdlib_symbols is \$(length(pkgnames_nonstdlib_symbols))\")",
        "println(\"The length of pkgnames_stdlib_symbols is \$(length(pkgnames_stdlib_symbols))\")",

        "println(\"pkgnames_symbols: (n = \$(length(pkgnames_symbols)))\")",
        "for i = 1:length(pkgnames_symbols)",
        "    println(\"\$(i). \$(pkgnames_symbols[i])\")",
        "end",

        "println(\"pkgnames_nonstdlib_symbols: (n = \$(length(pkgnames_nonstdlib_symbols)))\")",
        "for i = 1:length(pkgnames_nonstdlib_symbols)",
        "    println(\"\$(i). \$(pkgnames_nonstdlib_symbols[i])\")",
        "end",

        "println(\"pkgnames_stdlib_symbols: (n = \$(length(pkgnames_stdlib_symbols)))\")",
        "for i = 1:length(pkgnames_stdlib_symbols)",
        "    println(\"\$(i). \$(pkgnames_stdlib_symbols[i])\")",
        "end",

        "_simplecontainergenerator_containers_temp_depot_tmpdir_import_packagecompiler = mktempdir()",
        "atexit(() -> rm(_simplecontainergenerator_containers_temp_depot_tmpdir_import_packagecompiler; force = true, recursive = true))",
        "pushfirst!(Base.DEPOT_PATH, \"/opt/simplecontainergenerator_containers/packagecompiler_depot\")",
        "pushfirst!(Base.DEPOT_PATH, _simplecontainergenerator_containers_temp_depot_tmpdir_import_packagecompiler)",
        "import PackageCompiler",
        "popfirst!(Base.DEPOT_PATH)",
        "popfirst!(Base.DEPOT_PATH)",

        "make_sysimage = $(make_sysimage)",
        "if make_sysimage",
        "    PackageCompiler.create_sysimage(",
        "        pkgnames_symbols; ",
        "        cpu_target = \"$(julia_cpu_target)\", ",
        "        precompile_execution_file = \"/opt/simplecontainergenerator_containers/precompile_execution.jl\", ",
        "        project = \"/opt/simplecontainergenerator_containers/julia_project\", ",
        "        sysimage_path = \"/opt/simplecontainergenerator_containers/sysimage/SimpleContainerGeneratorSysimage.so\"",
        "        )",
        "end",
    ]
    content = join(lines, "\n") * "\n"
    return content
end


text:
using ModelingToolkit, OrdinaryDiffEq
using AnalogMTK
using AnalogMTK.EasyPlot #set, cons
import Printf: @sprintf
j=im

#==Circuit parameters
===============================================================================#
R = 1000
C = 10e-6 #10e-6, 10e-9
V = 1.0
𝑓 = 60 #Hz

#==Basic calculations
===============================================================================#
𝜔 = 2π*𝑓
ZC = 1/(j*𝜔*C)

println()
@info("Simulation parameters:")
@show V, 𝑓, 𝜔
@show R, C, abs(ZC)

𝜏 = R*C
@show 𝜏

#==Circuit definition
===============================================================================#
@named resistor = Resistor(R=R)
@named capacitor = Capacitor(C=C)
#@named source = VSource(V)
@named source = VSource(Sine(A=V, f=𝑓))
@named ground = Ground()

rc_eqs = [
          connect(source.p, resistor.p)
          connect(resistor.n, capacitor.p)
          connect(capacitor.n, source.n, ground.g)
         ]

@named rc_model = ODESystem(rc_eqs, AnalogMTK.t, systems=[resistor, capacitor, source, ground])

println()
@info("Initial ODESystem (model definition):"); flush(stdout); flush(stderr)
display(rc_model)

println()
@info("Computing simplified model (ODESystem)..."); flush(stdout); flush(stderr)
sys = structural_simplify(rc_model)
display(sys)
println("\tdone.")

#==Circuit simulation
===============================================================================#
println()
@info("Running simulation..."); flush(stdout); flush(stderr)
u0 = [
      capacitor.v => 0.0
      capacitor.p.i => 0.0
      resistor.v => 0.0
     ]
prob = ODEProblem(sys, u0, (0, 10.0))
sol = solve(prob, Trapezoid())
	#Trapezoid(), ImplicitEuler()
	#Rodas4(), RadauIIA5(), Tsit5()
@show length(sol.t)
println("\tdone.")

#==Plot results
===============================================================================#
println()
@info("Plotting results..."); flush(stdout); flush(stderr)
ylabels = ["Voltage [V]", "Current [A]"]
plot = plot_transient(sol.t, sol[source.p.v], id="Vsrc"; ylabels)
plot_transient!(plot, sol.t, sol[capacitor.p.v], id="Vout")
plot_transient!(plot, sol.t, sol[capacitor.p.i], id="Iout", strip=2)
displaygui(plot)
println("\tdone.")

:TEST_DONE


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

using Markdown
using InteractiveUtils

# ╔═╡ ae3abb60-13d1-11ec-0142-0fe32f29a9cf
begin
    using Pkg
    Pkg.activate(".")
    using Javis, Animations
end

# ╔═╡ c19adf8b-43af-4a55-b54f-de3b60855d83
begin
    function ground(args...)
        background("black") # canvas background
        sethue("white") # pen color
    end

    function object(p = O, color = "black")
        sethue(color)
        circle(p, 25, :fill)
        return p
    end

    function path!(points, pos, color)
        sethue(color)
        push!(points, pos) # add pos to points
        circle.(points, 2, :fill) # draws a circle for each point using broadcasting
    end

    function dancing_circles(c1, c2, start_pos = O)
        path_of_red = Point[]
        path_of_blue = Point[]

        red_ball = Object((args...) -> object(O, c1), start_pos + Point(100, 0))
        act!(red_ball, Action(anim_rotate_around(2π, start_pos)))
        blue_ball = Object((args...) -> object(O, c2), start_pos + Point(200, 80))
        act!(blue_ball, Action(anim_rotate_around(2π, 0.0, red_ball)))
        Object((args...) -> path!(path_of_red, pos(red_ball), c1))
        Object((args...) -> path!(path_of_blue, pos(blue_ball), c2))
    end
end

# ╔═╡ 51c3ad95-846b-420d-8be3-e82daf11ee5f
begin
    colored_planets = Video(500, 500)

    Background(1:70, ground)

    dancing_circles("green", "orange")

    render(colored_planets, pathname = "../output/colored_dancing_circles.gif")
end

# ╔═╡ 99beeaf4-1a48-4b7f-b48d-16b2705c1992
let

    moving_colored_planets = Video(500, 500)
    Background(1:140, ground)

    l1 = @JLayer 1:140 begin
        dancing_circles("green", "orange")
    end

    act!(l1, Action(1:1, anim_scale(0.4)))

    animation_point = Point(-75, -75)
    anim_back_and_forth = Animation(
        [0, 1 / 2, 1],
        [animation_point, -animation_point, animation_point],
        [sineio(), sineio()],
    )

    act!(l1, Action(1:140, anim_back_and_forth, translate()))

    render(
        moving_colored_planets,
        pathname = "../output/1_moving_colored_dancing_circles.gif",
    )
end

# ╔═╡ 3cc96cf6-1d12-4615-81af-f2529e0d7795
let

    finalvideo = Video(500, 500)

    Background(1:140, ground)

    colors = [
        ["red", "green"],
        ["orange", "blue"],
        ["yellow", "purple"],
        ["greenyellow", "darkgoldenrod1"],
    ]

    final_points = [Point(-150, -150), Point(150, -150), Point(150, 150), Point(-150, 150)]

    planets = map(colors) do c
        @JLayer 1:140 begin
            dancing_circles(c...)
        end
    end

    anim_back_and_forth = map(final_points) do point
        Animation([0.0, 1 / 2, 1.0], [O, point, O], [sineio(), sineio()])
    end

    map(zip(anim_back_and_forth, planets)) do (animation, pl)
        act!(pl, Action(1:1, anim_scale(0.4)))
        act!(pl, Action(1:140, animation, translate(), keep = false))
    end

    render(finalvideo; pathname = "../output/4_moving_colored_dancing_circles.gif")
end

# ╔═╡ Cell order:
# ╠═ae3abb60-13d1-11ec-0142-0fe32f29a9cf
# ╠═c19adf8b-43af-4a55-b54f-de3b60855d83
# ╠═51c3ad95-846b-420d-8be3-e82daf11ee5f
# ╠═99beeaf4-1a48-4b7f-b48d-16b2705c1992
# ╠═3cc96cf6-1d12-4615-81af-f2529e0d7795


text:
if basename(pwd()) == "aoc"
    cd("2016/23")
end
include("../12/12.jl")

expectedlength(line) = expectedlength(first(line))
function expectedlength(instr::AbstractString)
    if insorted(instr, ["dec", "inc", "tgl"])
        2
    elseif insorted(instr, ["add", "cpy", "jnz", "mul"])
        3
    else
        throw(ArgumentError("Unrecognized instruction: '$(instr)'"))
    end
end

function assembunny(filename::AbstractString; vm = AssembunnyVM(), pc = 1, debug = false, debugrun = false)
    lines = split.(readlines(filename))
    while checkbounds(Bool, lines, pc)
        line = lines[pc]
        pc += 1

        # step through each instruction
        if debug
            printstyled("pc: ", pc - 1, '\t', color = :blue)
            for sym in (:a, :b, :c, :d)
                printstyled(sym, ": "; color = :green)
                print(getfield(vm, sym), '\t')
            end
            printstyled(join(line, " "); color = :yellow)
            if debugrun
                println()
            else
                input = readline()
                if input == "stop"
                    return vm
                elseif input == "run"
                    debugrun = true
                end
            end
        end

        try
            if isempty(line) || length(line) != expectedlength(line)
                # do nothing

            elseif line[1] == "cpy"
                vm[line[3]] = line[2]

            elseif line[1] == "inc"
                if checkbounds(Bool, lines, pc) &&
                   checkbounds(Bool, lines, pc + 1) &&
                   lines[pc][1] == "dec" &&
                   lines[pc+1][1] == "jnz" &&
                   lines[pc][2] == lines[pc+1][2] &&
                   lines[pc+1][3] == "-2" &&
                   lines[pc][2] == lines[pc+1][2]

                    if checkbounds(Bool, lines, pc + 2) &&
                       checkbounds(Bool, lines, pc + 3) &&
                       lines[pc+2][1] == "dec" &&
                       lines[pc+3][1] == "jnz" &&
                       lines[pc+2][2] == lines[pc+3][2] &&
                       lines[pc+3][3] == "-5"

                        if debug
                            printstyled("OPTIMISED INSTR:\tmul ",
                                lines[pc+1][2], ' ',
                                lines[pc+3][2], ' ',
                                line[2], '\n';
                                color = :red)
                        end
                        vm[line[2]] += vm[lines[pc+1][2]] * vm[lines[pc+3][2]]
                        vm[lines[pc+1][2]] = 0
                        vm[lines[pc+3][2]] = 0
                        pc += 4

                    else
                        if debug
                            printstyled("OPTIMISED INSTR:\tadd ",
                                lines[pc+1][2], ' ', line[2], '\n'; color = :red)
                        end
                        vm[line[2]] += vm[lines[pc+1][2]]
                        vm[lines[pc+1][2]] = 0
                        pc += 2
                    end
                else
                    vm[line[2]] += 1
                end

            elseif line[1] == "dec"
                if lines[pc][1] == "dec" &&
                   lines[pc+1][1] == "jnz" &&
                   lines[pc+1][3] == "-2" &&
                   lines[pc][2] == lines[pc+1][2]

                    if debug
                        printstyled("OPTIMISED INSTR:\tSUB ", lines[pc+1][2], ' ', lines[pc][2], '\n'; color = :red)
                    end
                    vm[line[2]] -= vm[lines[pc][2]]
                    vm[lines[pc+1][2]] = 0
                    pc += 2
                else
                    vm[line[2]] -= 1
                end

            elseif line[1] == "jnz"
                if !iszero(vm[line[2]])
                    pc += vm[line[3]] - 1
                end

            elseif line[1] == "tgl"
                i = pc + vm[line[2]] - 1
                if checkbounds(Bool, lines, i)
                    instr = lines[i][1]
                    lines[i][1] =
                        if instr == "inc"
                            "dec"
                        elseif instr == "jnz"
                            "cpy"
                        elseif expectedlength(instr) == 2
                            "inc"
                        else
                            "jnz"
                        end
                end

            elseif line[1] == "add"
                vm[line[3]] += vm[line[2]]

            elseif line[1] == "sub"
                vm[line[3]] -= vm[line[2]]

            elseif line[1] == "mul"
                vm[line[3]] *= vm[line[2]]

            else
                throw(ArgumentError("Unrecognized instruction: '$(line[1])'"))

            end
        catch err
            println("Crashing on line ", pc - 1, ": ", join(line, ' '), '\t', (vm[x] for x in line[2:end])...)
            rethrow(err)
        end
    end
    if debug
        show(vm)
    end
    vm
end

function part1()
    @assert assembunny("example.txt").a == 3
    assembunny("input.txt"; vm = AssembunnyVM(7, 0, 0, 0)).a # 10365
end

function part2()
    assembunny("input.txt"; vm = AssembunnyVM(12, 0, 0, 0)).a
end

max_stars_repo_name: stellartux/advent-of-code

max_issues_repo_name: stellartux/advent-of-code

max_forks_repo_name: stellartux/advent-of-code


text:
const evtypes = String[ "not_existing",
                        "not_reported",
                        "anthropogenic_event",
                        "collapse",
                        "cavity_collapse",
                        "mine_collapse",
                        "building_collapse",
                        "explosion",
                        "accidental_explosion",
                        "chemical_explosion",
                        "controlled_explosion",
                        "experimental_explosion",
                        "industrial_explosion",
                        "mining_explosion",
                        "quarry_blast",
                        "road_cut",
                        "blasting_levee",
                        "nuclear_explosion",
                        "induced_or_triggered_event",
                        "rock_burst",
                        "reservoir_loading",
                        "fluid_injection",
                        "fluid_extraction",
                        "crash",
                        "plane_crash",
                        "train_crash",
                        "boat_crash",
                        "other_event",
                        "atmospheric_event",
                        "sonic_boom",
                        "sonic_blast",
                        "acoustic_noise",
                        "thunder",
                        "avalanche",
                        "snow_avalanche",
                        "debris_avalanche",
                        "hydroacoustic_event",
                        "ice_quake",
                        "slide",
                        "landslide",
                        "rockslide",
                        "meteorite",
                        "volcanic_eruption"]

"""
    randSeisSrc()

Generate a SeisSrc structure filled with random values.
"""
function randSeisSrc(; mw::Float32=0.0f0)
  R = SeisSrc()

#    src:
#     st: dur 0.0, rise 0.0, decay 0.0

  m0 = 1.0e-7*(10^(1.5*((mw == zero(Float32) ? 10^(rand(Float64)-1.0) : Float64(mw))+10.7)))

  setfield!(R, :id, string(rand(1:2^62)))                           # :id
  setfield!(R, :m0, m0)                                             # :m0
  setfield!(R, :gap, rand(Float64)*180.0)                           # :gap
  pop_rand_dict!(getfield(R, :misc), rand(4:24))                    # :misc
  setfield!(R, :mt, m0.*(rand(Float64, 6).-0.5))                    # :mt
  setfield!(R, :dm, m0.*(rand(Float64, 6).-0.5))                    # :dm
  setfield!(R, :npol, rand(Int64(6):Int64(120)))                    # :npol
  [note!(R, randstring(rand(16:256))) for i = 1:rand(1:6)]          # :notes
  setfield!(R, :pax,
    vcat(360.0.*(rand(Float64, 2, 3).-0.5),
          (rand(-1:2:1)*m0).*[rand() -1.0*rand() rand()]))          # :pax
  setfield!(R, :planes, 360.0.*(rand(Float64, 3, 2).-0.5))          # :planes
  setfield!(R, :src, "randSeisSrc")                                 # :src
  # setfield!(R, :src, join([join([randstring(12), " ",
  #                           randstring(20)]) for i=1:5], ","
  #                         ))                                        # :src
  setfield!(R, :st, SourceTime( randstring(2^rand(6:8)),            # :st
                                rand(Float64),
                                rand(Float64),
                                rand(Float64)))
  return R
end

"""
    randSeisHdr()

Generate a SeisHdr structure filled with random values.
"""
function randSeisHdr()
  H = SeisHdr()

  # moment magnitude and m0 in SI units
  mw = Float32(log10(10.0^rand(0.0:0.1:7.0)))
  m0 = 1.0e-7*(10^(1.5*(Float64(mw)+10.7)))

  # a random earthquake location
  dmin = 10.0^(rand()+0.5)
  dmax = dmin + exp(mw)*10.0^rand()

  loc = EQLoc(
           (rand(0.0:1.0:89.0) + rand())*-1.0^rand(1:2),
           (rand(0.0:1.0:179.0)+rand())*-1.0^rand(1:2),
           min(10.0*randexp(Float64), 660.0),
           rand()*4.0,
           rand()*4.0,
           rand()*8.0,
           rand(),
           rand(),
           rand(),
           10.0*10.0^rand(),
           dmin,
           dmax,
           4 + round(Int64, rand()*exp(mw)),
           rand(0x00:0x10:0xf0),
           rand_datum(),
           "hypocenter",
           "",
           rand(["HYPOELLIPSE", "HypoDD", "Velest", "centroid"]))

  # setfield!(loc, :sig, "1σ")
  # event type
  rtyp = rand()
  typ = rtyp > 0.3 ? "earthquake" : rand(evtypes)

  # Generate random header
  setfield!(H, :id, string(rand(1:2^62)))                           # :id
  setfield!(H, :int, (mw < 0.0f0 ? 0x00 : floor(UInt8, mw),
    randstring(rand(2:4))))                                         # :int
  setfield!(H, :loc, loc)                                           # :loc
  setfield!(H, :mag, EQMag( val   = mw,                             # :mag
                            scale = rand(["M","m"]) * randstring(2),
                            nst   = rand(3:100),
                            gap   = 360.0*rand(Float64),
                            src   = randstring(24)
                            )
            )
  # setfield!(H, :mech, FocalMech(src = join([join([randstring(12),     # :mech
  #                                     " ",
  #                                     randstring(20)]) for i=1:5], ","),
  #                               np  = 360.0.*(rand(Float64, 3, 2).-0.5),
  #                               pax = vcat(360.0.*(rand(Float64, 2, 2).-0.5),
  #                                     (rand(-1:2:1)*m0).*[rand() -1.0*rand()]),
  #                               m0  = m0,
  #                               mt  = m0.*(rand(Float64, 6).-0.5),
  #                               dm  = m0.*(rand(Float64, 6).-0.5)
  #                               )
  #         ) # obviously not pure DC, jajaja "random"
  pop_rand_dict!(H.misc, rand(4:24))                                # :misc
  [note!(H, randstring(rand(16:256))) for i = 1:rand(1:6)]          # :notes
  setfield!(H, :ot, now())                                          # :ot
  setfield!(H, :src, "randSeisHdr")                                 # :src
  setfield!(H, :typ, typ)                                           # :typ

  return H
end

"""
    randPhaseCat()

Generate a random seismic phase catalog suitable for testing EventChannel,
EventTraceData, and SeisEvent objects.
"""
function randPhaseCat(n::Int64=0)
  phase_list = String["P", "PKIKKIKP", "PKIKKIKS", "PKIKPPKIKP", "PKPPKP",
  "PKiKP", "PP", "PS", "PcP", "S", "SKIKKIKP", "SKIKKIKS", "SKIKSSKIKS", "SKKS",
  "SKS", "SKiKP", "SP", "SS", "ScS", "pP", "pPKiKP", "pS", "pSKS", "sP",
  "sPKiKP", "sS", "sSKS"]
  pol_list = Char['U', 'D', '-', '+', '_', ' ']

  phase = Array{String,1}(undef, n == 0 ? rand(3:18) : n)
  for j = 1:length(phase)
    phase[j] = rand(phase_list)
  end
  unique!(phase)

  P = PhaseCat()
  for j in phase
    P[j] = SeisPha(rand(), rand(), rand(), rand(), rand(), rand(), rand(), rand(), rand(pol_list), Char(rand(0x30:0x39)))
  end
  return P
end

"""
    randSeisEvent([, c=0.2, s=0.6])

Generate a SeisEvent structure filled with random header and channel data.
* 100*c is the percentage of :data channels _after the first_ with irregularly-sampled data (fs = 0.0)
* 100*s is the percentage of :data channels _after the first_ with guaranteed seismic data.
"""
function randSeisEvent(N::Int64; c::Float64=0.0, s::Float64=1.0, nx::Int64=0)
  V = SeisEvent(hdr=randSeisHdr(),
                source=randSeisSrc(),
                data=convert(EventTraceData, randSeisData(N, c=c, s=s, nx=nx)))

  for i = 1:V.data.n
    setindex!(getfield(getfield(V, :data), :pha), randPhaseCat(), i)
    setindex!(getfield(getfield(V, :data), :az),   180*(rand()-0.5), i)
    setindex!(getfield(getfield(V, :data), :baz),  180*(rand()-0.5), i)
    setindex!(getfield(getfield(V, :data), :dist), 180*rand(),       i)
  end
  return V
end
randSeisEvent(; c::Float64=0.0, s::Float64=1.0, nx::Int64=0) = randSeisEvent(rand(8:24), c=c, s=s, nx=nx)


text:
using Plots
plot(rand(10,10))

using SciMLBenchmarks
SciMLBenchmarks.bench_footer(WEAVE_ARGS[:folder],WEAVE_ARGS[:file])

max_stars_repo_name: Gregliest/SciMLBenchmarksOutput

max_issues_repo_name: Gregliest/SciMLBenchmarksOutput

max_forks_repo_name: Gregliest/SciMLBenchmarksOutput


text:
using Unitful
import Unitful: Frequency, Hz, kHz, 𝐓, s

samplerate = 48000Hz

abstract type AbstractOscillator end

function maxstable(samplerate::Frequency)
    return samplerate / 2 #TODO: this is not the correct formula
end

function checkstability(frequency::Frequency, samplerate::Frequency)
    if frequency >= maxstable(samplerate)
        error("Oscillator is unstable: $(frequency) is to high for $(samplerate) sample rate. The maximum stable frequency is $(maxstable(samplerate))")
    end
end

mutable struct SineOscillator{T<:AbstractFloat} <: AbstractOscillator
    previous::T
    sample::T
    cosine::T
end

function SineOscillator(frequency::Frequency{T}, amplitude::Complex{T}) where {T<:AbstractFloat}
    checkstability(frequency, samplerate)

    magnitude = abs(amplitude)
    phase = angle(amplitude)

    step = 2 * pi * frequency / samplerate
    previous = magnitude * sin(phase)
    sample = magnitude * sin(phase + step)
    cosine = cos(step)

    return SineOscillator{typeof(step)}(previous, sample, cosine)
end

function SineOscillator(frequency::Frequency{T}, amplitude::T) where {T<:AbstractFloat}
    return SineOscillator(frequency, complex(amplitude))
end

function oscillate!(osc::SineOscillator{T}) where {T<:AbstractFloat}
    nextsample = 2 * osc.sample * osc.cosine - osc.previous
    osc.previous, osc.sample = osc.sample, nextsample
    return nextsample
end

struct VectorOscillator{T<:AbstractFloat}
    oscs::Vector{SineOscillator{T}}
end

function oscillate!(osc::VectorOscillator{T}) where {T<:AbstractFloat}
    sum(oscillate!.(osc.oscs))
end


text:
using Multilane
using StatsBase
using MCTS
using JLD
using POMCP

dpws = DPWSolver(depth=20,
                 n_iterations=500,
                 exploration_constant=50.0,
                 k_state=4.0,
                 alpha_state=1/8,
                 rollout_solver=SimpleSolver())

pomcps = POMCPDPWSolver(
    eps=0.001,
    max_depth=20,
    c=50.0,
    tree_queries=2500,
    k_observation=4.0,
    alpha_observation=1/8,
    rollout_solver=SimpleSolver()
)

solvers = Dict{String, Any}(
    "dpw"=>dpws,
    "assume_normal"=>SingleBehaviorSolver(dpws, Multilane.NORMAL),
    "pomcp"=>MLPOMDPSolver(pomcps, BehaviorRootUpdaterStub(0.05))
)

curve = TestSet(lambda=[0.1, 1.0, 2.15, 4.64, 10., 21.5, 46.4], brake_threshold=2.5, N=500)
# curve = TestSet(lambda=[46.4], N=1)

tests = [
    TestSet(curve, solver_key="dpw", behaviors="9_agents", key="dpw"),
    TestSet(curve, solver_key="assume_normal", behaviors="9_agents", key="assume_normal"),
    TestSet(curve, solver_key="pomcp", behaviors="9_agents", key="pomcp")
]

objects = gen_initials(tests, generate_physical=true)

@show objects["param_table"]
objects["solvers"] = solvers

# files = sbatch_spawn(tests, objects,
#                      batch_size=50,
#                      time_per_batch="1:00:00",
#                      submit_command="sbatch",
#                      template_name="sherlock.sh")

results = evaluate(tests, objects, parallel=true)

filename = string("results_", Dates.format(Dates.now(),"u_d_HH_MM"), ".jld")
results["histories"] = nothing
save(filename, results)
println("results saved to $filename")


text:
# These tests don't do much except ensure that the package loads,
# and does something sensible if it does.
using Loess
using Test
using Random
using Statistics

Random.seed!(100)
xs = 10 .* rand(100)
ys = sin.(xs) .+ 0.5 * rand(100)

model = loess(xs, ys)

us = collect(minimum(xs):0.1:maximum(xs))
vs = predict(model, us)

@test minimum(vs) >= -1.1
@test maximum(vs) <= +1.1

x = [13.0,14.0,14.35,15.0,16.0]
y = [0.369486,  0.355579, 0.3545, 0.356952, 0.36883]
model = loess(x,y)
@test Loess.predict(model,x) ≈ y

let x = 1:10, y = sin.(1:10)
    model = loess(x, y)
    @test model.xs == reshape(collect(Float64, 1:10), 10, 1)
    @test model.ys == y
    pred = [1.02866, 0.798561, 0.533528, 0.253913, -0.0325918, -0.319578, -0.648763]
    @test predict(model, 1.0:0.5:4.0) ≈ pred atol=1e-5
end

@test_throws DimensionMismatch loess([1.0 2.0; 3.0 4.0], [1.0])


text:
# Parameters:
#	N::Int
#	- the problem size
#	K::Int
#	- the problem is homogenous at N-K
#	P::Array{Array{Array{Int, 1}, 1}, 1}
#	- output1 from partitions_bl()
#	ZFI::Array{Int, 1}
#	- output2 from partitions_bl()
#	WI::Array{Int, 2}
#	- output3 from partitions_bl()
# Return Values:
#	PT::Array{Array{Array{Int, 1}, 1}, 1} (Partition Tree)
#	- for each value, i, in PT[n][j], P[n][j] decomposes into P[n-1][i]
#	- length(PT[n]) = 0 for n = 1:(N - K)
#	- length(PT[n][p]) = 0 for p <= ZFI[n]
function partition_tree_bl(N::Int, K::Int, P::Array{Array{Array{Int, 1}, 1}, 1}, ZFI::Array{Int, 1}, WI::Array{Int, 2})
	PT = Array(Array{Array{Int, 1}, 1}, N)
	for n = 1:(N - K)
		PT[n] = Array(Array{Int, 1}, 0)
	end
	for n = (N - K + 1):N
		Pd = P[n - 1]
		Pn = P[n]
		PTn = Array(Array{Int, 1} , length(Pn))
		for p = 1:ZFI[n]
			PTn[p] = Array(Int, 0)
		end
		for p = (ZFI[n] + 1):length(Pn)
			PDA = pda(Pn[p])
			if length(PDA) == 1 || PDA[1][1] == PDA[2][1]
				lb = 1
				if PDA[1][1] != 1
					lb = WI[n - 1, PDA[1][1] - 1] + 1
				end
				PTn[p] = dia(PDA, Pd, lb)
			else
				lb1 = 1
				if PDA[1][1] != 1
					lb1 = WI[n - 1, PDA[1][1] - 1] + 1
				end
				lb2 = WI[n - 1, PDA[2][1] - 1]
				PTn[p] = dia(PDA, Pd, lb1, lb2)
			end
		end
		PT[n] = PTn
	end
	return PT
end


text:
"""
This evaluator is specialised for dealing with boolean multiplexer problems.
"""
type MultiplexerEvaluator <: SimpleEvaluator
  input_lines::Int
  select_lines::Int
  total_lines::Int
  permutations::Int
  inputs::Vector{Vector{Bool}}
  outputs::Vector{Bool}

  MultiplexerEvaluator(bits::Int) =
    build(new(2^bits, bits))
end

"""
Boolean multiplexer evaluator.

**Parameters:**

* `bits::Int`, the number of selection bits for this multiplexer.
"""
function multiplexer(bits::Int)
  e = MultiplexerEvaluator(bits)
  e.total_lines = e.input_lines + e.select_lines
  e.permutations = 2 ^ e.total_lines
  e.outputs = Array(Bool, e.permutations)
  e.inputs = Array(Vector{Bool}, e.permutations)

  for i in 1:e.permutations
    value = i
    divisor = e.permutations
    e.inputs[i] = Array(Bool, e.total_lines)

    for j in 1:e.total_lines
      divisor /= 2
      if value >= divisor
        e.inputs[i][j] = true
        value = divisor
      else
        e.inputs[i][j] = false
      end
    end

    choice = e.select_lines
    for (j, k) in enumerate(e.inputs[i][1:e.select_lines])
      if k
        choice += 2 ^ (j - 1)
      end
    end
    e.outputs[i] = e.inputs[i][choice]
  end

  return e
end

"""
function evaluate!(e::MultiplexerEvaluator, s::State, sc::FitnessScheme, c::Individual)
  outputs = Array(Bool, e.permutations)
  for i in 1:e.permutations
    inp = e.inputs[i]
    outputs[i] = execute(get(c.genome), inp[1], inp[2], inp[3],
      inp[4], inp[5], inp[6], inp[7], inp[8], inp[9], inp[10], inp[11])
  end

  hits = 0
  for i in 1:e.permutations
    if outputs[i] == e.outputs[i]
      hits += 1
    end
  end

  fitness(sc, hits)
end
"""


text:
using Flux
using Flux.Data: DataLoader
using Flux.Optimise: Optimiser, WeightDecay
using Flux: onehotbatch, onecold, glorot_normal, label_smoothing
using Flux.Losses: logitcrossentropy
using Statistics, Random
using Logging: with_logger
using ProgressMeter: @showprogress
import MLDatasets
using CUDA
using Formatting

using DeepUncertainty
include("utils.jl")

function LeNet5(args; imgsize = (28, 28, 1), nclasses = 10, device = gpu)
    out_conv_size = (imgsize[1] ÷ 4 - 3, imgsize[2] ÷ 4 - 3, 16)

    return Chain(
        SpectralConv((5, 5), imgsize[end] => 6, relu, device = device),
        MaxPool((2, 2)),
        SpectralConv((5, 5), 6 => 16, relu, device = device),
        MaxPool((2, 2)),
        flatten,
        SpectralDense(prod(out_conv_size), 120, relu, device = device),
        SpectralDense(120, 84, relu, device = device),
        SpectralDense(84, nclasses, device = device),
    )
end

loss(ŷ, y) = logitcrossentropy(ŷ, y)

function eval_loss_accuracy(args, loader, model, device)
    ntot = 0
    mean_l = 0
    mean_acc = 0
    mean_ece = 0
    mean_entropy = predictions = []
    targets = []
    for (x, y) in loader
        x, y = x |> device, y |> device
        model_preds = model(x)
        push!(predictions, cpu(model_preds))
        push!(targets, cpu(y))
        mean_l += loss(model_preds, y) * size(model_preds)[end]
        ntot += size(model_preds)[end]
    end
    # Calculate mean loss
    mean_l = mean_l / ntot |> round4
    predictions = Flux.batch(predictions)
    targets = Flux.batch(targets)
    pred_shape = size(predictions)
    target_shape = size(targets)
    predictions = reshape(predictions, (pred_shape[1], pred_shape[2] * pred_shape[3]))
    targets = reshape(targets, (target_shape[1], target_shape[2] * target_shape[3]))
    mean_acc = accuracy(predictions, targets) / ntot * 100 |> round4
    mean_entropy = round4(mean(calculate_entropy(predictions)))
    mean_ece = round4(expected_calibration_error(predictions, onecold(targets)))

    @info (format(
        "Mean Loss: {} Mean Accuracy: {} Mean ECE: {} Mean Entropy: {}",
        mean_l,
        mean_acc,
        mean_ece,
        mean_entropy,
    ))
    @info "==========================================================="
    return nothing
end

# arguments for the `train` function
Base.@kwdef mutable struct Args
    η = 3e-4             # learning rate
    λ = 0                # L2 regularizer param, implemented as weight decay
    batchsize = 256      # batch size
    epochs = 10          # number of epochs
    seed = 0             # set seed > 0 for reproducibility
    use_cuda = true      # if true use cuda (if available)
    infotime = 1      # report every `infotime` epochs
    checktime = 5        # Save the model every `checktime` epochs. Set to 0 for no checkpoints.
    dropout = 0.1
    sample_size = 10
end

function train(; kws...)
    args = Args(; kws...)
    args.seed > 0 && Random.seed!(args.seed)
    use_cuda = args.use_cuda && CUDA.functional()

    if use_cuda
        device = gpu
        @info "Training on GPU"
    else
        device = cpu
        @info "Training on CPU"
    end

    ## DATA
    train_loader, test_loader, ood_train_loader, ood_test_loader = get_data(args)
    @info "Dataset MNIST: $(train_loader.nobs) train and $(test_loader.nobs) test examples"

    ## MODEL AND OPTIMIZER
    model = LeNet5(args) |> device
    @info "LeNet5 model: $(num_params(model)) trainable params"

    ps = Flux.params(model)

    opt = ADAM(args.η)
    if args.λ > 0 # add weight decay, equivalent to L2 regularization
        opt = Optimiser(WeightDecay(args.λ), opt)
    end

    function report(epoch)
        @info "Test metrics"
        eval_loss_accuracy(args, test_loader, model, device)
        @info "Test OOD metrics"
        eval_loss_accuracy(args, ood_test_loader, model, device)
    end

    ## TRAINING
    @info "Start Training"
    report(0)
    for epoch = 1:args.epochs
        @showprogress for (x, y) in train_loader
            x, y = x |> device, y |> device
            gs = Flux.gradient(ps) do
                ŷ = model(x)
                total_loss = loss(ŷ, y)
                return total_loss
            end

            Flux.Optimise.update!(opt, ps, gs)
        end

        ## Printing and logging
        epoch % args.infotime == 0 && report(epoch)
    end
end

train()


text:
@testset "Frequency Adaptation" begin

    @testset "returns all indices once in the first block" begin
        for reps in 1:20
            n = rand(1:50)
            fa = BlackBoxOptim.FrequencyAdapter(n)
            block = Int[]
            for i in 1:n
                mi = BlackBoxOptim.next(fa)
                @test mi >= 1
                @test mi <= n
                push!(block, mi)
                BlackBoxOptim.update!(fa, mi, rand())
            end
            @test sort(block) == collect(1:n)
        end
    end

    @testset "increases the frequency of a method that has higher progress values" begin
        for reps in 1:20
            n = rand(2:50)
            fa = BlackBoxOptim.FrequencyAdapter(n)
            up1(mi) = BlackBoxOptim.update!(fa, mi, rand() + ((mi == 1) ? 0.4 : 0.0))
            counts = zeros(Int, n)
            uprep(rs) = begin
                for i in 1:rs
                    mi = BlackBoxOptim.next(fa)
                    counts[mi] += 1
                    up1(mi)
                end
            end
            uprep(n)
            @test counts == ones(Int, n)
            uprep(20*n)
            @test counts[1] > sum(counts[2:end])/(n-1)
        end
    end
end


text:
using ..TopOpt: GPU, @init_cuda
import ..TopOpt: whichdevice
import ..CUDASupport
using ..GPUUtils

@init_cuda

@define_cu(FilterMetadata, :cell_neighbouring_nodes, :cell_node_weights)
whichdevice(m::FilterMetadata) = whichdevice(m.cell_neighbouring_nodes)

@define_cu(SensFilter, :metadata, :nodal_grad, :last_grad, :cell_weights)
whichdevice(c::SensFilter) = whichdevice(c.nodal_grad)

function update_nodal_grad!(nodal_grad::CuVector, node_cells, args...)
    T = eltype(nodal_grad)
    allargs = (nodal_grad, node_cells.offsets, node_cells.values, args...)
    callkernel(dev, cheq_kernel1, allargs)
    CUDAdrv.synchronize(ctx)
    return
end

function cheq_kernel1(nodal_grad, node_cells_offsets, node_cells_values, cell_weights,
        cells, cellvolumes, black, white, varind, grad)
    T = eltype(nodal_grad)
    n = @thread_global_index()
    offset = @total_threads()
    while n <= length(nodal_grad)
        nodal_grad[n] = zero(T)
        cell_weights[n] = zero(T)
        r = node_cells_offsets[n]:node_cells_offsets[n+1]-1
        for i in r
            c = node_cells_values[i][1]
            if black[c] || white[c]
                continue
            end
            ind = varind[c]
            w = cellvolumes[c]
            cell_weights[n] += w
            nodal_grad[n] += w * grad[ind]
        end
        n += offset
    end
end

function normalize_grad!(nodal_grad::CuVector, cell_weights)
    T = eltype(nodal_grad)
    args = (nodal_grad, cell_weights)
    callkernel(dev, cheq_kernel2, args)
    CUDAdrv.synchronize(ctx)
    return
end
function cheq_kernel2(nodal_grad, cell_weights)
    T = eltype(nodal_grad)
    n = @thread_global_index()
    offset = @total_threads()
    while n <= length(nodal_grad)
        w = cell_weights[n]
        w = ifelse(w > 0, w, one(T))
        nodal_grad[n] /= w
        n += offset
    end
    return
end

function update_grad!(grad::CuVector, black, white, varind, cell_neighbouring_nodes, cell_node_weights, nodal_grad)
    T = eltype(grad)
    allargs = (grad, black, white, varind, cell_neighbouring_nodes.offsets, cell_neighbouring_nodes.values, cell_node_weights.values, nodal_grad)
    callkernel(dev, cheq_kernel3, allargs)
    CUDAdrv.synchronize(ctx)
    return
end
function cheq_kernel3(grad, black, white, varind, cell_neighbouring_nodes_offsets, cell_neighbouring_nodes_values, cell_node_weights_values, nodal_grad)
    T = eltype(nodal_grad)
    i = @thread_global_index()
    offset = @total_threads()
    while i <= length(black)
        if black[i] || white[i]
            continue
        end
        ind = varind[i]
        r = cell_neighbouring_nodes_offsets[ind]:cell_neighbouring_nodes_offsets[ind+1]-1
        length(r) == 0 && continue
        grad[ind] = zero(T)
        sum_weights = zero(T)
        for linear_ind in r
            w = cell_node_weights_values[linear_ind]
            sum_weights += w
            node_ind = cell_neighbouring_nodes_values[linear_ind]
            grad[ind] += nodal_grad[node_ind] * w
        end
        grad[ind] /= sum_weights
        i += offset
    end

    return
end


text:
using MuJoCo
using LinearAlgebra
using Printf
using Test

modelfile = joinpath(dirname(@__FILE__), "humanoid.xml")
#modelfile = dirname(@__FILE__)*"/swimmer.xml"

# write your own tests here
@info "Testing General Mujoco Environment"
include("test_env.jl")
@info "Testing Simple mjModel / mjData Access"
include("test_modeldata.jl")
#@info "Testing Finite Differencing Derivatives"
#include("test_deriv.jl")

#@info "Testing Name Access Method"
#include("test_names.jl")


text:
#****************************************************************************
#                     Molecular Dynamics Potentials (MDP)
#                            CESMIX-MIT Project
#
#  Contributing authors: Ngoc-Cuong Nguyen (cuongng@mit.edu, exapde@gmail.com)
#****************************************************************************

function genpotnum(filename, npm, sp0, sp1, sp2, sp3)

opufile = "opu" * filename;

tmp = "template <typename T> void " * opufile;
tmp = tmp * sp0;
tmp = tmp * "{\n";
for k = 1:npm
    if k == 1
        tmp = tmp * "\tif (potnum == " * string(k) * ")\n";
    else
        tmp = tmp * "\telse if (potnum == " * string(k) * ")\n";
    end
    tmp = tmp * "\t\t" * opufile * string(k) * sp3;
end
tmp = tmp * "}\n";
tmp = tmp * "template void " * opufile;
tmp = tmp * sp1;
tmp = tmp * "template void " * opufile;
tmp = tmp * sp2;

stropu = tmp;
strcpu = replace(stropu, "opu" => "cpu");
strgpu = replace(stropu, "opu" => "gpu");

return stropu, strcpu, strgpu

end

max_stars_repo_name: cesmix-mit/MDP

max_issues_repo_name: cesmix-mit/MDP

max_forks_repo_name: cesmix-mit/MDP


text:
using FuseApi
using Test

@testset "FuseApi.jl" begin
end


text:
n = [1,4,7,9,11,17,243,5247]

@test reinterpret(Int, decode(Base64, encode(Base64, reinterpret(UInt8, n)))) == n

for i in n
    @test decode(BCD, encode(BCD, i)) == i
    @test decode(BCD, encode(BCD, i, true), true) == i
end


text:
# Compute fast Walsh-Hadamard transforms (in natural, dyadic, or sequency order)
# using FFTW, by interpreting them as 2x2x2...x2x2 DFTs.  We follow Matlab's
# convention in which ifwht is the unnormalized transform and fwht has a 1/N
# normalization (as opposed to using a unitary normalization).
__precompile__(true)
module Hadamard
export fwht, ifwht, fwht_natural, ifwht_natural, fwht_dyadic, ifwht_dyadic, hadamard

using Base.FFTW
import Base.FFTW: set_timelimit, dims_howmany, unsafe_execute!, cFFTWPlan, r2rFFTWPlan, PlanPtr
import Base.DFT: normalization, complexfloat

using Compat

# A power-of-two dimension to be transformed is interpreted as a
# 2x2x2x....x2x2  multidimensional DFT.  This function transforms
# each individual power-of-two dimension n into the corresponding log2(n)
# DFT dimensions.  If bitreverse is true, then the output is in
# bit-reversed (transposed) order (which may not work in-place).
function hadamardize(dims::Array{Int,2}, bitreverse::Bool)
    ntot = 0
    for i = 1:size(dims,2)
        n = dims[1,i]
        if !ispow2(n)
            throw(ArgumentError("non power-of-two Hadamard-transform length"))
        end
        ntot += trailing_zeros(n)
    end
    hdims = Array(Int,3,ntot)
    j = 0
    for i = 1:size(dims,2)
        n = dims[1,i]
        is = dims[2,i]
        os = dims[3,i]
        log2n = trailing_zeros(n)
        krange = j+1:j+log2n
        for k in krange
            hdims[1,k] = 2
            hdims[2,k] = is
            hdims[3,k] = os
            is *= 2
            os *= 2
        end
        if bitreverse
            a = hdims[3,krange] # size = 1x3 in Julia 0.4, = 3 in 0.5
            hdims[3,krange] = flipdim(a, ndims(a))
        end
        j += log2n
    end
    return hdims
end

for (Tr,Tc,fftw,lib) in ((:Float64,:Complex128,"fftw",FFTW.libfftw),
                         (:Float32,:Complex64,"fftwf",FFTW.libfftwf))
    @eval function Plan_Hadamard{N}(X::StridedArray{$Tc,N}, Y::StridedArray{$Tc,N},
                                 region, flags::Unsigned, timelimit::Real,
                                 bitreverse::Bool)
        set_timelimit($Tr, timelimit)
        dims, howmany = dims_howmany(X, Y, [size(X)...], region)
        dims = hadamardize(dims, bitreverse)
        plan = ccall(($(string(fftw,"_plan_guru64_dft")),$lib),
                     PlanPtr,
                     (Int32, Ptr{Int}, Int32, Ptr{Int},
                      Ptr{$Tc}, Ptr{$Tc}, Int32, UInt32),
                     size(dims,2), dims, size(howmany,2), howmany,
                     X, Y, FFTW.FORWARD, flags)
        set_timelimit($Tr, NO_TIMELIMIT)
        if plan == C_NULL
            error("FFTW could not create plan") # shouldn't normally happen
        end
        return cFFTWPlan{$Tc,FFTW.FORWARD,X===Y,N}(plan, flags, region, X, Y)
    end

    @eval function Plan_Hadamard{N}(X::StridedArray{$Tr,N}, Y::StridedArray{$Tr,N},
                                 region, flags::Unsigned, timelimit::Real,
                                 bitreverse::Bool)
        set_timelimit($Tr, timelimit)
        dims, howmany = dims_howmany(X, Y, [size(X)...], region)
        dims = hadamardize(dims, bitreverse)
        kind = Array(Int32, size(dims,2))
        kind[:] = R2HC
        plan = ccall(($(string(fftw,"_plan_guru64_r2r")),$lib),
                     PlanPtr,
                     (Int32, Ptr{Int}, Int32, Ptr{Int},
                      Ptr{$Tr}, Ptr{$Tr}, Ptr{Int32}, UInt32),
                     size(dims,2), dims, size(howmany,2), howmany,
                     X, Y, kind, flags)
        set_timelimit($Tr, NO_TIMELIMIT)
        if plan == C_NULL
            error("FFTW could not create plan") # shouldn't normally happen
        end
        return r2rFFTWPlan{$Tr,(map(Int,kind)...),X===Y,N}(plan, flags, region, X, Y)
    end
end

############################################################################
# Define ifwht (inverse/unnormalized) transforms for various orderings

# Natural (Hadamard) ordering:

function ifwht_natural{T<:fftwNumber}(X::StridedArray{T}, region)
    Y = similar(X)
    p = Plan_Hadamard(X, Y, region, ESTIMATE, NO_TIMELIMIT, false)
    unsafe_execute!(p)
    return Y
end

function ifwht_natural{T<:Number}(X::StridedArray{T}, region)
    Y = float(X)
    p = Plan_Hadamard(Y, Y, region, ESTIMATE, NO_TIMELIMIT, false)
    unsafe_execute!(p)
    return Y
end

# Dyadic (Paley, bit-reversed) ordering:

function ifwht_dyadic{T<:fftwNumber}(X::StridedArray{T}, region)
    Y = similar(X)
    p = Plan_Hadamard(X, Y, region, ESTIMATE, NO_TIMELIMIT, true)
    unsafe_execute!(p)
    return Y
end

function ifwht_dyadic{T<:Number}(X::StridedArray{T}, region)
    return ifwht_dyadic(float(X), region)
end

############################################################################
# Sequency (Walsh) ordering:

# ifwht along a single dimension d of X
function ifwht{T<:fftwNumber}(X::Array{T}, region)
    Y = ifwht_dyadic(X, region)

    # Perform Gray-code permutation of Y (TODO: in-place?)
    if isempty(region)
        return Y
    elseif ndims(Y) == 1
        return [ Y[1 + ((i >> 1) $ i)] for i = 0:length(Y)-1 ]
    else
        sz = [size(Y)...]
        tmp = Array(T, maximum(sz[region])) # storage for out-of-place perm.
        for d in region
            # somewhat ugly loops to do 1d permutations along dimension d
            na = prod(sz[d+1:end])
            n = sz[d]
            nb = prod(sz[1:d-1])
            sa = nb * n
            for ia = 0:na-1
                for ib = 1:nb
                    i0 = ib + sa * ia
                    for i = 0:n-1
                        tmp[i+1] = Y[i0 + nb * ((i >> 1) $ i)]
                    end
                    for i = 0:n-1
                        Y[i0 + nb * i] = tmp[i+1]
                    end
                end
            end
        end
        return Y
    end
end

# handle 1d case of strided arrays (loops in multidim case are too annoying)
function ifwht{T<:fftwNumber}(X::StridedVector{T}, region)
    Y = ifwht_dyadic(X, region)

    # Perform Gray-code permutation of Y (TODO: in-place?)
    if isempty(region)
        return Y
    else
        return [ Y[1 + ((i >> 1) $ i)] for i = 0:length(Y)-1 ]
    end
end

# fallback for subarrays
function ifwht{T<:fftwNumber}(X::StridedArray{T}, region)
    return ifwht(copy(X), region)
end

# fallback for other types
function ifwht{T<:Number}(X::StridedArray{T}, region)
    return ifwht(float(X), region)
end

############################################################################
# Forward transforms (normalized by 1/N as in Matlab) and transforms
# without the region argument:

for f in (:ifwht_natural, :ifwht_dyadic, :ifwht)
    g = Symbol(string(f)[2:end])
    @eval begin
        $f(X) = $f(X, 1:ndims(X))
        $g(X) = scale!($f(X), normalization(X,1:ndims(X)))
        $g(X,r) = scale!($f(X,r), normalization(X,r))
    end
end

############################################################################
# Utilities to work with a precomputed cache of known Hadamard matrices
# of various sizes, produced by util/fetchhadamard.jl from Sloane's web page
# and stored as BitMatrices

function readcache(cachefile::AbstractString)
    B = BitMatrix[]
    open(cachefile, "r") do io
        while !eof(io)
            k = convert(Int, ntoh(read(io, Int64)))
            b = BitArray(k, k)
            # Hack: use internal binary data from BitArray for efficiency
            bits = read(io, eltype(b.chunks), length(b.chunks))
            for i = 1:length(bits)
                b.chunks[i] = ntoh(bits[i])
            end
            push!(B, b)
        end
    end
    return B
end
readcache() = readcache(joinpath(dirname(@__FILE__), "..", "src", "cache.dat"))

function printsigns{T<:Real}(io::IO, A::AbstractMatrix{T})
    m, n = size(A)
    println(io, m, "x", n, " sign matrix from ", typeof(A))
    for i = 1:m
        for j = 1:n
            print(io, A[i,j] > 0 ? "+" : "-")
        end
        println(io)
    end
end
printsigns(A) = printsigns(STDOUT, A)

function frombits(B::BitMatrix)
    A = convert(Matrix{Int8}, B)
    for i = 1:length(A)
        A[i] = A[i] == 0 ? -1 : 1
    end
    return A
end

const hadamards = BitMatrix[] # read lazily below

############################################################################
# Create Int8 order-n Hadamard matrix (in natural order), by factorizing
# n into a product of known Hadamard sizes (if possible)

const H2 = Int8[1 1; 1 -1]

function hadamard(n::Integer)
    n < 1 && throw(ArgumentError("size n=$n should be positive"))
    n0 = n
    H = reshape(Int8[1], 1,1)
    if !ispow2(n)
        isempty(hadamards) && append!(hadamards, readcache())
        for i = length(hadamards):-1:1
            k = size(hadamards[i], 1)
            if rem(n, k) == 0
                Hk = frombits(hadamards[i])
                while true
                    H = kron(H, Hk)
                    n = div(n, k)
                    rem(n, k) != 0 && break
                end
            end
        end
    end
    if !ispow2(n)
        n >>= trailing_zeros(n)
        throw(ArgumentError("unknown Hadamard factor $n of $n0"))
    end
    # Note: it would be faster to do a "power-by-squaring" like algorithm
    # here, where we repeatedly double the order via kron(H, H), but
    # it's not clear to me that we care much about performance here.
    while n > 1
        H = kron(H, H2)
        n >>= 1
    end
    return H
end

############################################################################

end # module


text:
immutable PolygonGeometry <: Gadfly.GeometryElement
    order::Int
    fill::Bool
    preserve_order::Bool

    function PolygonGeometry(; order::Int=0, fill::Bool=false,
                               preserve_order::Bool=false)
        return new(order, fill, preserve_order)
    end
end

const polygon = PolygonGeometry

function element_aesthetics(::PolygonGeometry)
    return [:x, :y, :color, :group]
end

function polygon_points(xs, ys, preserve_order)
    T = (@compat Tuple{eltype(xs), eltype(ys)})
    if preserve_order
        return T[(x, y)for (x, y) in zip(xs, ys)]
    else
        centroid_x, centroid_y = mean(xs), mean(ys)
        θ = atan2(xs - centroid_x, ys - centroid_y)
        perm = sortperm(θ)
        return T[(x, y)for (x, y) in zip(xs[perm], ys[perm])]
    end
end

# Render polygon geometry.
function render(geom::PolygonGeometry, theme::Gadfly.Theme,
                aes::Gadfly.Aesthetics)
    Gadfly.assert_aesthetics_defined("Geom.polygon", aes, :x, :y)

    default_aes = Gadfly.Aesthetics()
    default_aes.color = PooledDataArray(RGBA{Float32}[theme.default_color])
    aes = inherit(aes, default_aes)

    ctx = context(order=geom.order)
    T = (eltype(aes.x), eltype(aes.y))

    if aes.group != nothing
        XT, YT = eltype(aes.x), eltype(aes.y)
        xs = DefaultDict(Any, Vector{XT}, () -> XT[])
        ys = DefaultDict(Any, Vector{YT}, () -> YT[])
        for (x, y, c, g) in zip(aes.x, aes.y, cycle(aes.color), cycle(aes.group))
            push!(xs[(c,g)], x)
            push!(ys[(c,g)], y)
        end

        compose!(ctx,
            Compose.polygon([polygon_points(xs[(c,g)], ys[(c,g)], geom.preserve_order)
                             for (c,g) in keys(xs)]))
        cs = [c for (c,g) in keys(xs)]
        if geom.fill
            compose!(ctx, fill(cs),
                     stroke(map(theme.discrete_highlight_color, cs)))
        else
            compose!(ctx, fill(nothing), stroke(cs))
        end
    elseif length(aes.color) == 1 &&
            !(isa(aes.color, PooledDataArray) && length(levels(aes.color)) > 1)
        compose!(ctx, Compose.polygon(polygon_points(aes.x, aes.y, geom.preserve_order)))
        if geom.fill
            compose!(ctx, fill(aes.color[1]),
                     stroke(theme.discrete_highlight_color(aes.color[1])))
        else
            compose!(ctx, fill(nothing),
                     stroke(aes.color[1]))
        end
    else
        XT, YT = eltype(aes.x), eltype(aes.y)
        xs = DefaultDict(ColorValue, Vector{XT}, () -> XT[])
        ys = DefaultDict(ColorValue, Vector{YT}, () -> YT[])
        for (x, y, c) in zip(aes.x, aes.y, cycle(aes.color))
            push!(xs[c], x)
            push!(ys[c], y)
        end

        compose!(ctx,
            Compose.polygon([polygon_points(xs[c], ys[c], geom.preserve_order)
                             for c in keys(xs)]))
        cs = collect(keys(xs))
        if geom.fill
            compose!(ctx, fill(cs),
                     stroke(map(theme.discrete_highlight_color, cs)))
        else
            compose!(ctx, fill(nothing), stroke(cs))
        end
    end

    return compose!(ctx, linewidth(theme.line_width), svgclass("geometry"))
end


text:
# used as a notebook
using Plots, GLM, DataFrames
using Polynomials

xs = collect(5:5:50)
vs = [3.98, 6.3, 7.73, 8.66, 9.23, 9.46, 9.65, 9.86, 9.88, 9.89]

V₀ = 10
yp = log.(V₀ .- vs)

plot(xs, yp)

p = Polynomials.fit(xs, yp, 1)

tau = -1 / p[1]

b(n) = V₀ / (n * π) * (n % 2 == 0 ? 0 : 2)

T = 100
v_ft(t, N) = V₀ / 2 + sum([b(n) * sin(2π * n * t / T) for n ∈ 1:N])

v15 = [v_ft(t, 15) for t ∈ 5:5:40]
plot(collect(1:100), v15)

v100 = [v_ft(t, 100) for t ∈ 5:5:40]
plot!(collect(1:100), v100)

vs = vs[1:8]
yp = log.(v100 .- vs)
xs = collect(5:5:40)


text:
function GetScores(colorTree ::Array{Int,2})::Vector{Float64}  #::NamedTuple{(:alpha, :beta, :gamma)}
        numRows = size(colorTree, 1)
        numCols = size(colorTree, 2)
        if numCols<2 return NaN end
        global fnEvals+=1
        alpha = 0.0
        beta = 0.0
        gamma = 0.0
        for i=1:numRows
                #fromNode = colorTree[i,1]
                fromNode = i
                toNode = colorTree[i,1]
                if toNode>=fromNode
                        throw(ArgumentError("the colorTree is incorrectly formed as the toNode is larger than the fromNode"))
                end
                toColor = toNode==0 ? 7 : colorTree[toNode,2]
                fromColor = colorTree[fromNode,2]
                alpha += i * alphaScoreMatrix[fromColor, toColor]
                beta += i * betaScoreMatrix[fromColor, toColor]
                gamma += i * gammaScoreMatrix[fromColor, toColor]
        end
        return [alpha, beta, gamma]
        # return numRows .* [alpha, beta, gamma]
end

function GetCombinedScore(colorTree ::Array{Int,2})::Float64
        scores = GetScores(colorTree)
        return LinearAlgebra.norm(scores - target)
end

max_stars_repo_name: DesignEngrLab/ColorGraph

max_issues_repo_name: DesignEngrLab/ColorGraph

max_forks_repo_name: DesignEngrLab/ColorGraph


text:
@testset "GAN" begin
	using DrWatson
	@quickactivate
	using Test
	using GenerativeAD
	import StatsBase: fit!, predict
	using ConditionalDists, GenerativeModels, Distributions
	using EvalMetrics
	using Flux
	import GenerativeAD.Models: discriminate, generate

	# toy example
	M = 4
	N1 = 100
	N2 = 50
	X = vcat(randn(Float32, Int(M/2), N1), ones(Float32, Int(M/2), N1))
	Y = vcat(randn(Float32, Int(M/2), N2), ones(Float32, Int(M/2), N2))
	Z = vcat(ones(Float32, Int(M/2), N2), randn(Float32, Int(M/2), N2))
	data = ((X,),(Y,zeros(size(Y,2))))
	parameters = (zdim=8, hdim=32, lr=0.001, batchsize=8, activation="swish", nlayers=3, idim=M)
	function basic_convergence_test(model, gloss, dloss; kwargs...)
		history, iterations, model = fit!(model, data, gloss, dloss; patience = 100, max_iter=5000,
			max_train_time=600, kwargs...)
		scores = map(x->GenerativeAD.Models.discriminate(model, x), (X,Y,Z))
		@test mean(scores[1]) > mean(scores[3])
		@test mean(scores[2]) > mean(scores[3])
		model
	end
	# vanilla GAN
	model = GenerativeAD.Models.gan_constructor(;parameters...)
	gloss = GenerativeAD.Models.gloss
	dloss = GenerativeAD.Models.dloss
	model = basic_convergence_test(model, gloss, dloss; parameters...)
	@test size(GenerativeAD.Models.generate(model, 10)) == (M,10)
	# feature-matching GAN
	model = GenerativeAD.Models.gan_constructor(;parameters...)
	alpha = 1f0
	fmloss(args...) = alpha*GenerativeAD.Models.gloss(args...) + GenerativeAD.Models.fmloss(args...)
	dloss = GenerativeAD.Models.dloss
	model = basic_convergence_test(model, fmloss, dloss; parameters...)
	# wasserstein gan
	model = GenerativeAD.Models.gan_constructor(;last_linear=true,parameters...)
	dps = Flux.params(model.discriminator)
	wgloss(model,x) = - mean(discriminate(model, generate(model, size(x,ndims(x)))))
	wdloss(model,x) = - mean(discriminate(model,x)) + mean(discriminate(model, generate(model, size(x,ndims(x)))))
	#model = basic_convergence_test(model, wgloss, wdloss; weight_clip=0.01,
	#	discriminator_advantage=5, stop_threshold=-Inf, max_iter=1500, parameters...)
end


text:
using Distributions
using Statistics

println("analysis here")

max_stars_repo_name: banyan-team/banyan-julia

max_issues_repo_name: banyan-team/banyan-julia

max_forks_repo_name: banyan-team/banyan-julia


text:
slurp(s) = readcsv(IOBuffer(s))

conv(s)= colon(map(x->parse(Int,x),match(r"^(-?\d+)-(-?\d+)$", s).captures)...)

expand(s) = mapreduce(x -> isa(x,Number)? Int(x) : conv(x), vcat, slurp(s))

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
function vtkWriteQuad(filename,X,Connectivity,c)

precision = '2';

fid=fopen(filename, 'w');
fprintf(fid, '# vtk DataFile Version 2.0 \n');
fprintf(fid, 'Unstructured Grid Example \n');
fprintf(fid, 'ASCII \n');
fprintf(fid, 'DATASET UNSTRUCTURED_GRID \n');
fprintf(fid, 'POINTS %d float \n',size(X,2));

spec = ['%0.', precision, 'f '];
fprintf(fid, spec, X);
fprintf(fid, '\n');
fprintf(fid, 'CELLS %d %d \n',size(Connectivity,2),5*size(Connectivity,2));
fprintf(fid,'4 %d %d %d %d \n',Connectivity-1);
fprintf(fid, 'CELL_TYPES %d \n',size(Connectivity,2));
fprintf(fid,'%d  \n',9*ones(size(Connectivity,2),1));
fprintf(fid, 'Cell_DATA %d \n',size(Connectivity,2));
fprintf(fid, 'SCALARS scalars float 1 \n');
fprintf(fid, 'LOOKUP_TABLE default \n');
fprintf(fid,'%d %d %d %d \n',c);
fclose(fid);
end


text:
using Test
using Distributions

include("../kernels/kernel.jl")
include("../computation/finitedifference.jl")
include("../computation/tdist.jl")
include("../transforms/transforms.jl")
include("../datastructs.jl")

@testset "kernel.jl" begin
    @test BoxCox()(3, 1) ≈ 2.0 atol = 1e-9
    @test BoxCox()(inverse(BoxCox(), 2 ,3), 3) ≈ 2 atol = 1e-9
    #@test invBoxCox(boxCox(2, 3), 3) ≈ 2.0 atol = 1e-15
    #@test invBoxCox(boxCox(2, 1), 1) ≈ 2.0 atol = 1e-15
    #@test boxCox(invBoxCox(2, 3), 3) ≈ 2.0 atol = 1e-15
    #@test boxCox(invBoxCox(2, 1), 1) ≈ 2.0 atol = 1e-15
    #@test boxCox(invBoxCox([1 2 3], 1), 1) ≈ [1 2 3] atol = 1e-15
    #@test boxCox(invBoxCox([1; 2], 1), 1) ≈ [1;2] atol = 1e-15
    #@test boxCox(invBoxCox(2., -0.45), -0.45) ≈ 2.0 atol = 1e-15
    #@test invBoxCox(boxCox(2., -2), -2) ≈ 2.0 atol = 1e-15
end

@testset "BoxCox" begin
    bc = BoxCox()
    g = (x, λ) -> bc(x, λ)
    dg = (x, λ) -> partialx(bc, x, λ)
    dg2 = (x, λ) -> partialx(bc, x, λ)
    g_fixed = x-> g(x, 1.5)
    dg_fixed = x-> dg(x, 1.5)
    dg2_fixed = x-> dg2(x, 1.5)
    (_, _, _, pol) = checkDerivative(dg_fixed, dg2_fixed, 0.5, nothing, 8, 16, 10)
    @test coeffs(pol)[end] ≈ 2 atol = 1e-1
    (_, _, _, pol) = checkDerivative(g_fixed, dg_fixed, 0.5, nothing, 8, 16, 10)
    @test coeffs(pol)[end] ≈ 2 atol = 1e-1
end

@testset "datastructs.jl" begin
    x = newExtensibleTrainingData(5, 2, 3)
    @test x.d == 2
    @test x.p == 3
    @test x.n == 0
    o = [1.0 2.0; 3.0 4.0]
    p = [1.0 2.0 3.0; 4.0 5.0 6.0]
    z = [1.0, 2.0]
    update!(x, o, p, z)
    @test x.x[1:2, :] == o
    @test x.Fx[1:2, :] == p
    @test x.y[1:2] == z
    @test x.n == 2
    o2 = [1.0 2.0; 3.0 4.0; 4.0 5.0]
    p2 = [1.0 2.0 3.0; 4.0 5.0 6.0; 1.0 2.0 3.0]
    z2 = [1.0;2.0; 3.0]
    update!(x, o2, p2, z2)
    @test x.x[1:5, :] == vcat(o, o2)
    @test x.Fx[1:5, :] == vcat(p, p2)
    @test x.y[1:5] == vcat(z, z2)
    @test x.n == 5
    o3 = [1.0 2.0]
    p3 = [1.0 2.0 3.0]
    z3 = [3.0]
   # @test_broken BoundsError update!(x, o3, p3, z3)
end

#@testset "kernel.jl" begin
#a = (tau, x) -> exp.(-x*tau/2)
#da = (tau, x) -> -tau/2 * exp.(-x*tau/2)
#gg = x -> a(2, x)
#dgg =  x -> da(2, x)
#(h, A, plt1, poly, fis1) = checkDerivative(gg, dgg, 1.0, nothing, 1, 2, 10)
#println(poly)
#display(plt1)

#dsdf = RBF()
#dsdf = Gaussian()
#ff = y -> dsdf(2, y)
#dff =  y -> partial_θ(dsdf, 2, y)
#(h, A, plt1, poly, fis2, vals, debugvals) = checkDerivative(ff, dff, 1.0, nothing, 1, 2, 10)
#println(poly)
#display(plt1)
#Plots.plot(vals[:, 1], seriestype = :scatter)
#for i =2:10
#    println("here")
#    Plots.plot!(vals[:, i], seriestype = :scatter)
#end
#Plots.plot(fis1')
#Plots.plot(fis2')
#plt(ff, .01, 10)
#plt!(dff, .01, 10)
#plt!(gg, .01, 10)
#plt!(dgg, .01, 10)
#end

max_stars_repo_name: dbindel/btg

max_issues_repo_name: dbindel/btg

max_forks_repo_name: dbindel/btg


text:
# TODO: clean up ugly routine
function translate(cursor::CLVarDecl)
    child_cursors = children(cursor)
    cursor_sym = Symbol(spelling(cursor))
    if isempty(child_cursors)
        return MetaExpr(cursor_sym, cursor)
    elseif length(child_cursors) == 1
        var_cursor = first(child_cursors)
        var_meta = translate(var_cursor)
        var_meta.expr == Expr(:null) && return MetaExpr(cursor_sym, cursor)
        var_meta.expr = Expr(:(=), cursor_sym, var_meta.expr)
        return var_meta
    else
        var_meta = translate(first(child_cursors))
        list_meta = translate(last(child_cursors))
        if Meta.isexpr(var_meta.expr, :call) && kind(last(child_cursors)) != CXCursor_CallExpr
            push!(var_meta.expr.args, list_meta.expr.args...)
            var_meta.expr = Expr(:(=), cursor_sym, var_meta.expr)
        else
            var_meta.expr = Expr(:(=), cursor_sym, list_meta.expr)
        end
        return var_meta
    end
end

function translate(cursor::CLTypeRef)
    origin = getref(cursor)
    if kind(origin) == CXCursor_StructDecl
        return MetaExpr(Expr(:call, Symbol(spelling(origin))))
    else
        # TODO: check this out
        @warn "unexpected subroutine reached."
        return MetaExpr(Expr(:null), cursor)
    end
end

translate(cursor::CLParmDecl) = MetaExpr(Symbol(spelling(cursor)), cursor)

function translate(cursor::CLFunctionDecl)
    child_cursors = children(cursor)
    signature = Expr(:call, Symbol(spelling(cursor)))
    body = Expr(:null)
    for c in child_cursors
        meta = translate(c)
        if kind(c) == CXCursor_ParmDecl
            push!(signature.args, meta.expr)
        else
            body = meta.expr
        end
    end
    return MetaExpr(Expr(:function, signature, body))
end


text:
# This is part of Quiz.jl. License is MIT.

include("types/abstract.jl")
include("types/union.jl")
include("types/qElem.jl")


text:
"""
A Julia package for the exploration of continuous
and discrete dynamical systems.
"""
module DynamicalSystems

"""
    DynamicalSystem
Abstract type representing a dynamical system. Has the following concrete sub-types:
* `DiscreteDS1D`
* `DiscreteDS`
* `ContinuousDS`
"""
abstract type DynamicalSystem end

export DynamicalSystem, Systems

include("discrete.jl")
include("continuous.jl")
include("lyapunovs.jl")
include("famous_systems.jl")

end # module


text:
module Promises

module 💛
include("./notebook.jl")
end

const Promise = 💛.Promise
const all = 💛.promise_all
const any = 💛.promise_any
const race = 💛.promise_race
const resolve = 💛.Promises.resolve
const reject = 💛.Promises.reject
const await_settled = 💛.await_settled
const Resolved = 💛.Resolved
const Rejected = 💛.Rejected
const PromiseSettledResult = 💛.PromiseSettledResult
const var"@await" = 💛.var"@await"
const var"@async_promise" = 💛.var"@async_promise"

export Promise,
    @await,
    @async_promise

end


text:
module BaryRational
export FHInterp, bary, chebpts, chebwts, floater_weights, lagrange_weights, aaa, prz, differ, Dcheb
using LinearAlgebra
using Statistics
using SparseArrays

abstract type BRInterp <: Function end

struct FHInterp{T <: AbstractArray} <: BRInterp
    x::T
    f::T
    w::T
    order::Int
end

include("weights.jl")
include("bary.jl")
include("aaa.jl")
include("derivatives.jl")

function FHInterp(x::Vector{T}, f::Vector{T}; order::Int=0, grid=false) where {T}
    if grid
        FHInterp(collect(x), f, floater_weights(length(x), order), order)
    else
        FHInterp(collect(x), f, floater_weights(x, order), order)
    end
end

(y::FHInterp)(z) = bary(z, y.f, y.x, y.w)

end