text:
abstract type Solver end

abstract type ForwardBackwardSolver <: Solver end

export Solver, ForwardBackwardSolver

################################################################################
export PG, FPG

"""
`PG(;kwargs...)`

Creates an object `PG` containing the options of the Proximal Gradient solvers:

  * `gamma`, stepsize (default: unspecified, determined automatically)
  * `maxit`, maximum number of iteration (default: `10000`)
  * `tol`, halting tolerance on the fixed-point residual (default: `1e-4`)
  * `adaptive`, adaptively adjust `gamma` (default: `false` if `gamma` is provided)
  * `fast`, enables accelerated method (default: `false`)
  * `verbose`, verbosity level (default: `1`)
  * `verbose_freq`, verbosity frequency for `verbose = 1` (default: `100`)

"""
struct PG <: ForwardBackwardSolver
    kwargs::Iterators.Pairs
    function PG(; kwargs...)
        new(kwargs)
    end
end

"""
`FPG(;kwargs...)`

Same as `PG`, creates the options of the Fast Proximal Gradient solver.

"""
function FPG(; kwargs...)
    return PG(; kwargs..., fast=true)
end

function build_iterator(x, solver::PG; kwargs...)
    x, ProximalAlgorithms.FBSIterator(~x; solver.kwargs..., kwargs...)
end

################################################################################
export ZeroFPR

"""
`ZeroFPR(;kwargs...)`

Creates an object `ZeroFPR` containing the options of the ZeroFPR solver:

  * `gamma`, stepsize (default: unspecified, determined automatically)
  * `maxit`, maximum number of iteration (default: `10000`)
  * `tol`, halting tolerance on the fixed-point residual (default: `1e-4`)
  * `adaptive`, adaptively adjust `gamma` (default: `false` if `gamma` is provided)
  * `fast`, enables accelerated method (default: `false`)
  * `verbose`, verbosity level (default: `1`)
  * `verbose_freq`, verbosity frequency for `verbose = 1` (default: `100`)
  * `memory`, memory of the `LBFGS` operator (default: `10` )

"""
struct ZeroFPR <: ForwardBackwardSolver
    kwargs::Iterators.Pairs
    function ZeroFPR(; kwargs...)
        new(kwargs)
    end
end

function build_iterator(x, solver::ZeroFPR; kwargs...)
    x, ProximalAlgorithms.ZeroFPRIterator(~x; solver.kwargs..., kwargs...)
end

################################################################################
export PANOC

"""
`ZeroFPR(;kwargs...)`

Creates an object `PANOC` containing the options of the PANOC solver:

  * `gamma`, stepsize (default: unspecified, determined automatically)
  * `maxit`, maximum number of iteration (default: `10000`)
  * `tol`, halting tolerance on the fixed-point residual (default: `1e-4`)
  * `adaptive`, adaptively adjust `gamma` (default: `false` if `gamma` is provided)
  * `fast`, enables accelerated method (default: `false`)
  * `verbose`, verbosity level (default: `1`)
  * `verbose_freq`, verbosity frequency for `verbose = 1` (default: `100`)
  * `memory`, memory of the `LBFGS` operator (default: `10` )

"""
struct PANOC <: ForwardBackwardSolver
    kwargs::Iterators.Pairs
    function PANOC(; kwargs...)
        new(kwargs)
    end
end

function build_iterator(x, solver::PANOC; kwargs...)
    x, ProximalAlgorithms.PANOCIterator(~x; solver.kwargs..., kwargs...)
end

default_solver = PANOC

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


text:
using Hecke

bounddisc = ARGS[1]
gtype = ARGS[2]

file = "conductors_" * bounddisc * "_" * gtype;

bounddisc = fmpz(eval(parse(bounddisc)))
gtype = convert(Vector{Int}, eval(parse(gtype)))

sprint_formatted(fmt, args...) = @eval @sprintf($fmt, $(args...))

Qx, x = PolynomialRing(QQ, "x")
K, a = NumberField(x - 1, "a")
O = maximal_order(K)

n=prod(gtype)
expo=lcm(gtype)

print("Computing conductors ... ")

conductors = sort!(Hecke.conductorsQQ(O, n, bounddisc))

println("DONE")

if length(ARGS) == 3
  file = ARGS[3]
else
  file = file * "_" * string(length(conductors))
end

if isfile(file)
  error("File $file already exists")
end

open(file, "w") do f
  print("Writing conductors to $file ... ")
  writedlm(f, conductors)
  println("DONE")
end

@assert countlines(file) == length(conductors)

ext: jl


text:
import Distributions: Normal
import UncertainData: UncertainValue, UncertainValueDataset, resample

# Vectors of uncertain values
uvals_x = [UncertainValue(Normal, rand(Normal(0, 5)), abs(rand(Normal(0, 3)))) for i = 1:100]
uvals_y = [UncertainValue(Normal, rand(Normal(0, 5)), abs(rand(Normal(0, 3)))) for i = 1:100];

# UncertainValueDataset
UVX = UncertainValueDataset(uvals_x)
UVY = UncertainValueDataset(uvals_y)

# UncertainIndexDataset
UVX_idx = UncertainIndexDataset(uvals_x)
UVY_idx = UncertainIndexDataset(uvals_y)

# Real-valued vectors
x = resample.(uvals_x);
y = resample.(uvals_y);

################################################################
# Integration with UncertainData.jl
################################################################
ctest = CrossMappingTest(n_reps = 10)

# We should get a vector of cross map values
@test causality(x, y, ctest) isa Vector{T} where T
@test causality(uvals_x, uvals_y, ctest) isa Vector{T} where T
@test causality(x, uvals_y, ctest) isa Vector{T} where T
@test causality(uvals_x, y, ctest) isa Vector{T} where T
@test causality(UVX, UVY, ctest) isa Vector{T} where T
@test causality(x, UVY, ctest) isa Vector{T} where T
@test causality(UVX, y, ctest) isa Vector{T} where T

# n_reps = 10, so we should get 10 cross map values
@test causality(x, y, ctest) |> length == 10
@test causality(uvals_x, uvals_y, ctest) |> length == 10
@test causality(x, uvals_y, ctest) |> length == 10
@test causality(uvals_x, y, ctest) |> length == 10
@test causality(UVX, UVY, ctest) |> length == 10
@test causality(x, UVY, ctest) |> length == 10
@test causality(UVX, y, ctest) |> length == 10

################################################################
# Integration with UncertainData.jl (with sampling constraints)
################################################################
onevar_constraints = ConstrainedResampling(TruncateStd(1))
twovar_constraints = ConstrainedResampling(TruncateStd(2), TruncateStd(1))

test_cm = CrossMappingTest()
@test causality(x, y, test_cm, twovar_constraints) isa Vector{T} where T<:Real
@test causality(uvals_x, uvals_y, test_cm, twovar_constraints) isa Vector{T} where T<:Real
@test causality(x, uvals_y, test_cm, onevar_constraints) isa Vector{T} where T<:Real
@test causality(uvals_x, y, test_cm, onevar_constraints) isa Vector{T} where T<:Real
@test causality(UVX, UVY, test_cm, twovar_constraints) isa Vector{T} where T<:Real
@test causality(uvals_x, UVY, test_cm, twovar_constraints) isa Vector{T} where T<:Real
@test causality(UVX, uvals_y, test_cm, twovar_constraints) isa Vector{T} where T<:Real
@test causality(x, UVY, test_cm, onevar_constraints) isa Vector{T} where T<:Real
@test causality(UVX, y, test_cm, onevar_constraints) isa Vector{T} where T<:Real


text:
using Param
using Base.Test

@paramdef(LVPar,
    r = 1.5,
    K = 10.0,
    a = 1.8,
    e = 1.0,
    m = 0.3
)

@withparam LVPar par function ydot(t, y, par)
    ydot = zeros(2)
    ydot[1] = r*y[1]*(1.0 - y[1]/K) - a*y[1]*y[2]
    ydot[2] = e*a*y[1]*y[2]
    return ydot
end

p1 = LVPar(r = 0.8, a = 1.1)
ydot(0, [0.5, 0.5], p1)

max_stars_repo_name: gabrielgellner/Param

max_issues_repo_name: gabrielgellner/Param

max_forks_repo_name: gabrielgellner/Param


text:
using Registrator: Registrator
using Registrator.WebUI: @gf
using GitForge: GitForge, get_user
using GitForge.GitHub: GitHub, GitHubAPI, NoToken, Token
using HTTP: HTTP
using Sockets: Sockets
using Distributed

const UI = Registrator.WebUI

empty!(UI.CONFIG)
merge!(UI.CONFIG, Dict(
    "ip" => "localhost",
    "port" => 4000,
    "registry_url" => "https://github.com/JuliaRegistries/General",
    "server_url" => "http://localhost:4000",
    "github" => Dict{String, Any}(
        # We need a token to avoid rate limits on Travis.
        "token" => get(ENV, "GITHUB_API_TOKEN", ""),
        "client_id" => "",
        "client_secret" => "",
    ),
    "gitlab" => Dict{String, Any}(
        "token" => "",
        "client_id" => "",
        "client_secret" => "",
    ),
))

const backup = deepcopy(UI.CONFIG)

function restoreconfig!()
    empty!(UI.CONFIG)
    merge!(UI.CONFIG, deepcopy(backup))
end

@testset "Web UI" begin
    @testset "Provider initialization" begin
        UI.init_providers()
        @test length(UI.PROVIDERS) == 2
        @test Set(collect(keys(UI.PROVIDERS))) == Set(["github", "gitlab"])
        empty!(UI.PROVIDERS)

        delete!(UI.CONFIG, "github")
        UI.CONFIG["gitlab"]["disable_rate_limits"] = true
        UI.init_providers()
        @test collect(keys(UI.PROVIDERS)) == ["gitlab"]
        @test !GitForge.has_rate_limits(UI.PROVIDERS["gitlab"].client, identity)
        empty!(UI.PROVIDERS)
        restoreconfig!()

        mktemp() do path, io
            extra_provider = """
                PROVIDERS["myprovider"] = Provider(;
                    name="MyProvider",
                    client=GitHubAPI(),
                    client_id="abc",
                    client_secret="abc",
                    auth_url="abc",
                    token_url="abc",
                    scope="public_repo",
                )
                """
            print(io, extra_provider)
            close(io)

            UI.CONFIG["extra_providers"] = path
            UI.init_providers()
            @test haskey(UI.PROVIDERS, "myprovider")
            empty!(UI.PROVIDERS)
            restoreconfig!()
        end
    end

    # Patch the GitHub API client to avoid needing a real API key.
    t = isempty(UI.CONFIG["github"]["token"]) ? NoToken() : Token(UI.CONFIG["github"]["token"])
    UI.init_providers()
    UI.PROVIDERS["github"] = UI.Provider(;
        name="GitHub",
        client=GitHubAPI(; token=t),
        client_id=UI.CONFIG["github"]["client_id"],
        client_secret=UI.CONFIG["github"]["client_secret"],
        auth_url="https://github.com/login/oauth/authorize",
        token_url="https://github.com/login/oauth/access_token",
        scope="public_repo",
    )

    @testset "Registry initialization" begin
        UI.init_registry()
        reg = UI.REGISTRY[]
        @test reg.url == reg.clone == UI.CONFIG["registry_url"]
        @test reg.repo isa GitHub.Repo

        UI.CONFIG["registry_clone_url"] = "git@github.com:JuliaRegistries/General.git"
        UI.init_registry()
        reg = UI.REGISTRY[]
        @test reg.url == UI.CONFIG["registry_url"]
        @test reg.clone == UI.CONFIG["registry_clone_url"]
        restoreconfig!()
    end

    # Start the server.
    # TODO: Stop it when this test set is done.
    task = @async UI.start_server(Sockets.localhost, 4000)
    @info "Waiting for server to start..."
    sleep(10)    # Wait for server to be up

    @testset "404s" begin
        for r in values(UI.ROUTES)
            resp = HTTP.get(UI.CONFIG["server_url"] * r * "/foo"; status_exception=false)
            @test resp.status == 404
            @test occursin("Page not found", String(resp.body))
        end
    end

    @testset "Route: /" begin
        # The response should contain the registry URL and authentication links.
        resp = HTTP.get(UI.CONFIG["server_url"]; status_exception=false)
        @test resp.status == 200
        s = String(resp.body)
        @test occursin(UI.CONFIG["registry_url"], s)
        @test occursin("GitHub", s)
        @test occursin("GitLab", s)

        # When a provider is disabled, its authentication link should not appear.
        delete!(UI.PROVIDERS, "gitlab")
        resp = HTTP.get(UI.CONFIG["server_url"]; status_exception=false)
        @test resp.status == 200
        s = String(resp.body)
        @test occursin("GitHub", s)
        @test !occursin("GitLab", s)
    end

    @testset "Route: /select" begin
        resp = HTTP.get(UI.CONFIG["server_url"] * UI.ROUTES[:SELECT]; status_exception=false)
        @test resp.status == 200
        @test occursin("package to register", String(resp.body))
    end

    @testset "Route: /register (validation)" begin
        resp = HTTP.get(UI.CONFIG["server_url"] * UI.ROUTES[:REGISTER]; status_exception=false)
        @test resp.status == 405
        @test occursin("Method not allowed", String(resp.body))

        resp = HTTP.post(UI.CONFIG["server_url"] * UI.ROUTES[:REGISTER]; status_exception=false)
        @test resp.status == 400
        @test occursin("Missing or invalid state cookie", String(resp.body))

        # Pretend we've gone through authentication.
        state = "foo"
        client = UI.PROVIDERS["github"].client
        user = @gf get_user(client, "octocat")
        UI.USERS[state] = UI.User(user, client)

        url = UI.CONFIG["server_url"] * UI.ROUTES[:REGISTER]
        cookies = Dict("state" => state)

        body = "package=&ref=master"
        resp = HTTP.post(url; body=body, cookies=cookies, status_exception=false)
        @test resp.status == 400
        @test occursin("Package URL was not provided", String(resp.body))

        body = "package=foo&ref=master"
        resp = HTTP.post(url; body=body, cookies=cookies, status_exception=false)
        @test resp.status == 400
        @test occursin("Package URL is invalid", String(resp.body))

        body = "package=https://github.com/foo/bar&ref="
        resp = HTTP.post(url; body=body, cookies=cookies, status_exception=false)
        @test resp.status == 400
        @test occursin("Branch was not provided", String(resp.body))

        body = "package=https://github.com/JuliaLang/NotARealRepo&ref=master"
        resp = HTTP.post(url; body=body, cookies=cookies, status_exception=false)
        @test resp.status == 400
        @test occursin("Repository was not found", String(resp.body))

        body = "package=http://github.com/JuliaLang/Example.jl&ref=master"
        resp = HTTP.post(url; body=body, cookies=cookies, status_exception=false)
        @test resp.status == 400
        @test occursin("Unauthorized to release this package", String(resp.body))

        body = "package=git@github.com:JuliaLang/Example.jl.git&ref=master"
        resp = HTTP.post(url; body=body, cookies=cookies, status_exception=false)
        @test resp.status == 400
        @test occursin("Unauthorized to release this package", String(resp.body))
    end

end


text:
mutable struct Stmt <: DBInterface.Statement
    con::Connection
    handle::duckdb_prepared_statement
    sql::AbstractString

    function Stmt(con::Connection, sql::AbstractString)
        handle = Ref{duckdb_prepared_statement}()
        if duckdb_prepare(con.handle, sql, handle) != DuckDBSuccess
            error_message = unsafe_string(duckdb_prepare_error(handle))
            duckdb_destroy_prepare(handle)
            throw(QueryException(error_message))
        end
        stmt = new(con, handle[], sql)
        finalizer(_close_stmt, stmt)
        return stmt
    end

    function Stmt(db::DB, sql::AbstractString)
        return Stmt(db.main_connection, sql)
    end
end

function _close_stmt(stmt::Stmt)
    if stmt.handle != C_NULL
        duckdb_destroy_prepare(stmt.handle)
    end
    stmt.handle = C_NULL
    return
end

DBInterface.getconnection(stmt::Stmt) = stmt.con

duckdb_bind_internal(stmt::Stmt, i::Integer, val::AbstractFloat) = duckdb_bind_double(stmt.handle, i, Float64(val));
duckdb_bind_internal(stmt::Stmt, i::Integer, val::Bool) = duckdb_bind_boolean(stmt.handle, i, val);
duckdb_bind_internal(stmt::Stmt, i::Integer, val::Int8) = duckdb_bind_int8(stmt.handle, i, val);
duckdb_bind_internal(stmt::Stmt, i::Integer, val::Int16) = duckdb_bind_int16(stmt.handle, i, val);
duckdb_bind_internal(stmt::Stmt, i::Integer, val::Int32) = duckdb_bind_int32(stmt.handle, i, val);
duckdb_bind_internal(stmt::Stmt, i::Integer, val::Int64) = duckdb_bind_int64(stmt.handle, i, val);
duckdb_bind_internal(stmt::Stmt, i::Integer, val::UInt8) = duckdb_bind_uint8(stmt.handle, i, val);
duckdb_bind_internal(stmt::Stmt, i::Integer, val::UInt16) = duckdb_bind_uint16(stmt.handle, i, val);
duckdb_bind_internal(stmt::Stmt, i::Integer, val::UInt32) = duckdb_bind_uint32(stmt.handle, i, val);
duckdb_bind_internal(stmt::Stmt, i::Integer, val::UInt64) = duckdb_bind_uint64(stmt.handle, i, val);
duckdb_bind_internal(stmt::Stmt, i::Integer, val::Float32) = duckdb_bind_float(stmt.handle, i, val);
duckdb_bind_internal(stmt::Stmt, i::Integer, val::Float64) = duckdb_bind_double(stmt.handle, i, val);
duckdb_bind_internal(stmt::Stmt, i::Integer, val::Date) = duckdb_bind_date(stmt.handle, i, value_to_duckdb(val));
duckdb_bind_internal(stmt::Stmt, i::Integer, val::Time) = duckdb_bind_time(stmt.handle, i, value_to_duckdb(val));
duckdb_bind_internal(stmt::Stmt, i::Integer, val::DateTime) =
    duckdb_bind_timestamp(stmt.handle, i, value_to_duckdb(val));
duckdb_bind_internal(stmt::Stmt, i::Integer, val::Missing) = duckdb_bind_null(stmt.handle, i);
duckdb_bind_internal(stmt::Stmt, i::Integer, val::Nothing) = duckdb_bind_null(stmt.handle, i);
duckdb_bind_internal(stmt::Stmt, i::Integer, val::AbstractString) = duckdb_bind_varchar(stmt.handle, i, val);
duckdb_bind_internal(stmt::Stmt, i::Integer, val::Vector{UInt8}) = duckdb_bind_blob(stmt.handle, i, val, sizeof(val));
duckdb_bind_internal(stmt::Stmt, i::Integer, val::WeakRefString{UInt8}) =
    duckdb_bind_varchar(stmt.handle, i, val.ptr, val.len);

function duckdb_bind_internal(stmt::Stmt, i::Integer, val::Any)
    println(val)
    throw(NotImplementedException("unsupported type for bind"))
end

function bind_parameters(stmt::Stmt, params::DBInterface.StatementParams)
    i = 1
    for param in params
        if duckdb_bind_internal(stmt, i, param) != DuckDBSuccess
            throw(QueryException("Failed to bind parameter"))
        end
        i += 1
    end
end

max_stars_repo_name: lokax/duckdb

max_issues_repo_name: lokax/duckdb

max_forks_repo_name: lokax/duckdb


text:
"Split a subspace into two"
split(ss::UInt64) = hash((ss, 0)), hash((ss, 1))

max_stars_repo_name: zenna/expect

max_issues_repo_name: zenna/expect

max_forks_repo_name: zenna/expect


text:
"""
    QuantizedArray{Q<:AbstractQuantization,U<:Unsigned,D<:PreMetric,T,N} <: AbstractArray{T,N}

A quantized array. It represents a 'quantized' representation of an original
array, equivalent to a lossy compressed version.

# Fields
  * `quantizer::ArrayQuantizer{Q,U,D,T,N}` the quantized of the array
  * `data::Matrix{U}` the actual compressed representation of the quantized array
"""
struct QuantizedArray{Q,U<:Unsigned,D<:Distances.PreMetric,T,N} <: AbstractArray{T,N}
    quantizer::ArrayQuantizer{Q,U,D,T,N}
    data::Matrix{U}
end

# Aliases
const QuantizedVector{Q,U,D,T} = QuantizedArray{Q,U,D,T,1}
const QuantizedMatrix{Q,U,D,T} = QuantizedArray{Q,U,D,T,2}
const OrthogonalQuantizedArray{U,D,T,N} = QuantizedArray{OrthogonalQuantization,U,D,T,N}
const OrthogonalQuantizedMatrix{U,D,T,N} = QuantizedArray{OrthogonalQuantization,U,D,T,2}
const AdditiveQuantizedArray{U,D,T,N} = QuantizedArray{AdditiveQuantization,U,D,T,N}
const AdditiveQuantizedMatrix{U,D,T,N} = QuantizedArray{AdditiveQuantization,U,D,T,2}

# Main outer constructor
function QuantizedArray(aa::AbstractArray{T,N};
                        k::Int=DEFAULT_K,
                        m::Int=DEFAULT_M,
                        method::Symbol=DEFAULT_METHOD,
                        distance::Distances.PreMetric=DEFAULT_DISTANCE,
                        kwargs...) where {T,N}
    @assert N <=2 "Array quantization is supported only for Vectors and Matrices"
    @assert k >= 1 "`k` has to be larger or equal to 1"
    @assert m >= 1 "`m` has to be larger or equal to 1"
    method != :rvq && @assert rem(nvars(aa), m) == 0 "`m` has to divide exactly $(nvars(aa))"
    method != :sample && @assert T<:AbstractFloat "Input array has to have AbstractFloat elements"
    aq = build_quantizer(aa, k=k, m=m, method=method, distance=distance; kwargs...)
    data = quantize_data(aq, aa)
    return QuantizedArray(aq, data)
end

"""
    quantize(aq, aa)

Quantize an array `aa` using an array quantizer `aq`.
"""
function quantize(aq::ArrayQuantizer{Q,U,D,T,N}, aa::AbstractArray{T,N}) where {Q,U,D,T,N}
    new_aq = ArrayQuantizer(aq.quantization, size(aa), codebooks(aq), aq.k, aq.distance, aq.rot)
    data = quantize_data(new_aq, aa)
    return QuantizedArray(new_aq, data)
end

"""
    quantize(aa; [;kwargs])
Quantize an array `aa`.

# Arguments
  * `aa::AbstractArray` input array to be quantized

# Keyword arguments
  * `k::Int` the number of vector prototypes in each codebook
  * `m::Int` the number of codebooks
  * `method::Symbol` the algorithm to be employed for codebook
generation; possible values are `:sample` (default), `:pq` for
classical k-means clustering codebooks and `:opq` for 'cartesian'
k-means clustering codebooks
  * `distance::PreMetric` the distance to be used in the
codebook generation methods and data encoding

Other codebook generation algorithm specific keyword arguments
such as `maxiter::Int` can be specified as well.
"""
quantize(aa::AbstractArray{T,N}; kwargs...) where {T,N} = QuantizedArray(aa; kwargs...)

# eltype, size, length
Base.eltype(qa::QuantizedArray{Q,U,D,T,N}) where {Q,U,D,T,N} = T
Base.size(qa::QuantizedArray) = qa.quantizer.dims
Base.IndexStyle(::Type{<:QuantizedArray}) = IndexLinear()

"""
    quantizer(qa)

Access the quantizer field of a quantized array `qa`.
"""
quantizer(qa::QuantizedArray) = qa.quantizer

"""
    nvars(aa)

Returns the number of variables of an array `aa`. For vectors,
the value is always `1`, for matrices it is the number of
rows in the matrix.
"""
nvars(av::AbstractVector) = 1
nvars(am::AbstractMatrix) = size(am, 1)

# Indexing interface: getindex, setindex!
@inline function Base.getindex(qa::OrthogonalQuantizedArray, i::Int)
    cbooks = codebooks(quantizer(qa))
    m = length(cbooks)
    col, row, cidx = _indices(nvars(qa), m, i)
    qkey = getindex(qa.data, cidx, col)
    return cbooks[cidx][qkey][row]
end

@inline function Base.getindex(qa::OrthogonalQuantizedMatrix, i::Int, j::Int)
    cbooks = codebooks(quantizer(qa))
    m = length(cbooks)
    _, row, cidx = _indices(nvars(qa), m, i)
    qkey = getindex(qa.data, cidx, j)
    return cbooks[cidx][qkey][row]
end

@inline function Base.getindex(qa::AdditiveQuantizedArray, i::Int)
    cbooks = codebooks(quantizer(qa))
    m = length(cbooks)
    col, row = divrem(i[1]-1, nvars(qa)) .+ 1
    return sum(cbooks[cidx][getindex(qa.data, cidx, col)][row] for cidx in 1:m)
end

@inline function Base.getindex(qa::AdditiveQuantizedMatrix, i::Int, j::Int)
    cbooks = codebooks(quantizer(qa))
    m = length(cbooks)
    return sum(cbooks[cidx][getindex(qa.data, cidx, j)][i] for cidx in 1:m)
end

function _indices(n::Int, m::Int, I::Vararg{Int,N}) where {N}
    step = Int(n/m)
    col, r = divrem(I[1]-1, n) .+ 1          # column, row in uncompressed data
    cidx, row = divrem(r-1, step) .+ 1       # codebook index, codebook vector row
    return col, row, cidx
end

# setindex! not supported
Base.setindex!(qa::QuantizedArray, i::Int) =
    throw(ErrorException("setindex! not supported on QuantizedArrays"))

Base.setindex!(qa::QuantizedArray, I::Vararg{Int, N}) where {N}=
    throw(ErrorException("setindex! not supported on QuantizedArrays"))


text:
function gamess_to_julia_key(gamess_key)
  gamess_to_julia_mapping = Dict([ 1 => 1,
                                   2 => 2,
                                   3 => 3,
                                   4 => 4,
                                   5 => 5,
                                   6 => 6,
                                   7 => 7,
                                   8 => 8,
                                   9 => 9,
                                   10 => 10,
                                   11 => 13,
                                   12 => 15,
                                   13 => 11,
                                   14 => 12,
                                   15 => 14,
                                   16 => 16,
                                   17 => 17,
                                   18 => 18,
                                   19 => 19,
                                   20 => 20
                                 ])

  ij_match = eachmatch(r"[0-9]+", gamess_key)
  ij = collect(ij_match)
  i = parse(Int,ij[1].match)
  j = parse(Int,ij[2].match)

  #i_new = gamess_to_julia_mapping[i]
  #j_new = gamess_to_julia_mapping[j]

  i_new = i
  j_new = j

  if i_new <= j_new
    tmp = i_new
    i_new = j_new
    j_new = tmp
  end

  return "$i_new,$j_new"
  #return "$i,$j,$k,$l"
end

function julia_to_gamess_key(julia_key)
  julia_to_gamess_mapping = Dict([ 1 => 1,
                                   2 => 2,
                                   3 => 3,
                                   4 => 4,
                                   5 => 5,
                                   6 => 6,
                                   7 => 7,
                                   8 => 8,
                                   9 => 9,
                                   10 => 10,
                                   13 => 11,
                                   15 => 12,
                                   11 => 13,
                                   12 => 14,
                                   14 => 15,
                                   16 => 16,
                                   17 => 17,
                                   18 => 18,
                                   19 => 19,
                                   20 => 20
                                 ])

  ij_match = eachmatch(r"[0-9]+", julia_key)
  ij = collect(ij_match)
  i = parse(Int,ij[1].match)
  j = parse(Int,ij[2].match)

  #i_new = julia_to_gamess_mapping[i]
  #j_new = julia_to_gamess_mapping[j]

  i_new = i
  j_new = j

  if i_new <= j_new
    tmp = i_new
    i_new = j_new
    j_new = tmp
  end

  return "$i_new,$j_new"
end

function get_matrix(input_file_string::String, matrix_string::String)
  #== read in input file ==#
  input_file::Array{String,1} = []
  open(input_file_string) do file
    input_file = readlines(file)
  end

  #== extract matrixs from input file ==#
  starting::Int64 = 1
  ending::Int64 = 1

  matrix = Vector{String}([])
  for line in input_file
    if occursin("$matrix_string(",line)
      push!(matrix,line)
    end
  end

  #== set up matrix dictionary ==#
  matrix_dict = Dict([])
  num_duplicates = 0
  for element in matrix
    #==get ij numbers ==#
    ij_match = eachmatch(r"[0-9]+", element)
    ij = collect(ij_match)
    i = ij[1].match
    j = ij[2].match

    #== get actual matrix ==#
    matrix_match = match(r"0","0")
    if typeof(match(r"(-?[\d]{1,}\.[\d]{1,}[e()E][+()-][\d]{1,})",element)) != Nothing
      matrix_match = match(r"(-?[\d]{1,}\.[\d]{1,}[e()E][+()-][\d]{1,})",element)
    elseif typeof(match(r"(-\.[\d]{1,})", element)) != Nothing
      matrix_match = match(r"(-\.[\d]{1,})",element)
    else
      matrix_match = match(r"(-?[\d]\.[\d]{1,})",element)
    end

    matrix_value = parse(Float64,matrix_match.match)
    if haskey(matrix_dict,"$i,$j")
      println("Duplicate key at: "*"$i,$j")
      num_duplicates += 1
    else
      matrix_dict["$i,$j"] = matrix_value
    end
  end

  println("Number of duplicate keys: ", num_duplicates)

  return matrix_dict
end

julia_dict = get_matrix(ARGS[1], ARGS[3])
gamess_dict = get_matrix(ARGS[2], ARGS[3])

println(length(keys(julia_dict)))
println(length(keys(gamess_dict)))

for gamess_key in keys(gamess_dict)
  julia_key = gamess_to_julia_key(gamess_key)
  has_julia_key = haskey(julia_dict,julia_key)
  if !has_julia_key
    println("Missing Julia key:", julia_key)
  else
    same = gamess_dict[gamess_key] < 1.0E-10 && julia_dict[julia_key] < 1.0E-10 ? true :
      isapprox(gamess_dict[gamess_key], julia_dict[julia_key], atol = 1e-9)
    if !same
      println("Not same!: $gamess_key => $julia_key: ", gamess_dict[gamess_key], ", ", julia_dict[julia_key])
    end
  end
end

for julia_key in keys(julia_dict)
  gamess_key = julia_to_gamess_key(julia_key)
  has_gamess_key = haskey(gamess_dict,gamess_key)
  if !has_gamess_key
    println("Missing GAMESS key:", gamess_key)
  else
    same = gamess_dict[gamess_key] < 1.0E-10 && julia_dict[julia_key] < 1.0E-10 ? true :
      isapprox(gamess_dict[gamess_key], julia_dict[julia_key], atol = 1e-9)
    if !same
      println("Not same!: $julia_key => $gamess_key: ", julia_dict[julia_key], ", ", gamess_dict[gamess_key])
    end
  end
end


text:
abstract type LenseFlowOp{I<:ODESolver,t₀,t₁,Φ} <: FlowOpWithAdjoint{I,t₀,t₁} end

# `L = LenseFlow(ϕ)` just creates a wrapper holding ϕ. Then when you do `L*f` or
# `cache(L,f)` we create a CachedLenseFlow object which holds all the
# precomputed quantities and preallocated memory needed to do the lense.

struct LenseFlow{I<:ODESolver,t₀,t₁,Φ<:Field{<:Any,<:S0}} <: LenseFlowOp{I,t₀,t₁,Φ}
    ϕ::Φ
end

struct CachedLenseFlow{N,t₀,t₁,Φ<:Field,ŁΦ<:Field,ÐΦ<:Field,ŁF<:Field,ÐF<:Field,T} <: LenseFlowOp{RK4Solver{N},t₀,t₁,Φ}

    # save ϕ to know when to trigger recaching
    ϕ :: Ref{Any}

    # p and M⁻¹ quantities precomputed at every time step
    p   :: Dict{Float16,FieldOrOpVector{Diagonal{T,ŁΦ}}}
    M⁻¹ :: Dict{Float16,FieldOrOpMatrix{Diagonal{T,ŁΦ}}}

    # f type memory
    memŁf  :: ŁF
    memÐf  :: ÐF
    memŁvf :: FieldVector{ŁF}
    memÐvf :: FieldVector{ÐF}

    # ϕ type memory
    memŁϕ  :: ŁΦ
    memÐϕ  :: ÐΦ
    memŁvϕ :: FieldVector{ŁΦ}
    memÐvϕ :: FieldVector{ÐΦ}

end

### constructors
LenseFlow(ϕ,n=7) = LenseFlow{RK4Solver{n}}(ϕ)
LenseFlow{I}(ϕ) where {I<:ODESolver} = LenseFlow{I,0,1}(ϕ)
LenseFlow{I,t₀,t₁}(ϕ) where {I,t₀,t₁} = LenseFlow{I,float(t₀),float(t₁),typeof(ϕ)}(ϕ)

### printing
show(io::IO, ::L) where {I,t₀,t₁,Φ,L<:LenseFlow{I,t₀,t₁,Φ}} = print(io, "$(L.name.name){$t₀→$t₁, $I}(ϕ::$Φ)")
show(io::IO, ::L) where {N,t₀,t₁,Φ,ŁF,L<:CachedLenseFlow{N,t₀,t₁,Φ,<:Any,<:Any,ŁF}} = print(io, "$(L.name.name){$t₀→$t₁, $(RK4Solver{N})}(ϕ::$Φ, Łf::$ŁF)")
string(::Type{RK4Solver{N}}) where {N} = "$N-step RK4"

# convenience for getting the actual ϕ map
getϕ(L::LenseFlow) = L.ϕ
getϕ(L::CachedLenseFlow) = L.ϕ[]

### caching
τ(t) = Float16(t)
cache(L::LenseFlow, f) = cache!(alloc_cache(L,f),L,f)
cache(cL::CachedLenseFlow, f) = cL
cache!(cL::CachedLenseFlow{N,t₀,t₁}, ϕ) where {N,t₀,t₁} =
    (cL.ϕ[]===ϕ) ? cL : cache!(cL,LenseFlow{RK4Solver{N},t₀,t₁}(batch_promote!(cL.memÐϕ,Ð(ϕ))),cL.memŁf)
(cL::CachedLenseFlow)(ϕ) = cache!(cL,ϕ)
function cache!(cL::CachedLenseFlow{N,t₀,t₁}, L::LenseFlow{RK4Solver{N},t₀,t₁}, f) where {N,t₀,t₁}
    ts = range(t₀,t₁,length=2N+1)
    ∇ϕ,Hϕ = map(Ł, gradhess(L.ϕ))
    T = eltype(L.ϕ)
    for (t,τ) in zip(ts,τ.(ts))
        @! cL.M⁻¹[τ] = pinv(Diagonal.(I + T(t)*Hϕ))
        @! cL.p[τ] = cL.M⁻¹[τ]' * Diagonal.(∇ϕ)
    end
    cL.ϕ[] = L.ϕ
    cL
end
function alloc_cache(L::LenseFlow{RK4Solver{N},t₀,t₁}, f) where {N,t₀,t₁}
    ts = range(t₀,t₁,length=2N+1)
    p, M⁻¹ = Dict(), Dict()
    Łf,Ðf = Ł(f),  Ð(f)
    Łϕ,Ðϕ = Ł(L.ϕ),Ð(L.ϕ)
    for (t,τ) in zip(ts,τ.(ts))
        M⁻¹[τ] = Diagonal.(similar.(@SMatrix[Łϕ Łϕ; Łϕ Łϕ]))
        p[τ]   = Diagonal.(similar.(@SVector[Łϕ,Łϕ]))
    end
    CachedLenseFlow{N,t₀,t₁,typeof(L.ϕ),typeof(Łϕ),typeof(Ðϕ),typeof(Łf),typeof(Ðf),eltype(Łϕ)}(
        Ref{Any}(L.ϕ), p, M⁻¹,
        similar(Łf), similar(Ðf), similar.(@SVector[Łf,Łf]), similar.(@SVector[Ðf,Ðf]),
        similar(Łϕ), similar(Ðϕ), similar.(@SVector[Łϕ,Łϕ]), similar.(@SVector[Ðϕ,Ðϕ]),
    )
end

# the way these velocities work is that they unpack the preallocated fields
# stored in L.mem* into variables with more meaningful names, which are then
# used in a bunch of in-place (eg mul!, Ł!, etc...) functions. note the use of
# the @! macro, which just rewrites @! x = f(y) to x = f!(x,y) for easier
# reading.

function velocity(L::LenseFlowOp{<:RK4Solver}, f₀::Field)
    function v!(v::Field, t::Real, f::Field)
        Ðf, Ð∇f, Ł∇f = L.memÐf, L.memÐvf,  L.memŁvf
        p = L.p[τ(t)]

        @! Ðf  = Ð(f)
        @! Ð∇f = ∇ᵢ * Ðf
        @! Ł∇f = Ł(Ð∇f)
        @! v   = p' * Ł∇f
    end
    return (v!, batch_promote!(L.memŁf,Ł(f₀)))
end

function velocityᴴ(L::LenseFlowOp{<:RK4Solver}, f₀::Field)
    function v!(v::Field, t::Real, f::Field)
        Łf, Łf_p, Ð_Łf_p = L.memŁf, L.memŁvf, L.memÐvf
        p = L.p[τ(t)]

        @! Łf = Ł(f)
        @! Łf_p = p * Łf
        @! Ð_Łf_p = Ð(Łf_p)
        @! v = -∇ᵢ' * Ð_Łf_p
    end
    return (v!, batch_promote!(L.memÐf,Ð(f₀)))
end

function negδvelocityᴴ(L::LenseFlowOp{<:RK4Solver}, (f₀, δf₀)::FieldTuple)

    function v!((df_dt, dδf_dt, dδϕ_dt)::FieldTuple, t::Real, (f, δf, δϕ)::FieldTuple)

        p   = L.p[τ(t)]
        M⁻¹ = L.M⁻¹[τ(t)]

        # dδf/dt
        Łδf, Łδf_p, Ð_Łδf_p = L.memŁf, L.memŁvf, L.memÐvf
        @! Łδf     = Ł(δf)
        @! Łδf_p   = p * Łδf
        @! Ð_Łδf_p = Ð(Łδf_p)
        @! dδf_dt  = -∇ᵢ' * Ð_Łδf_p

        # df/dt
        Ðf, Ð∇f, Ł∇f = L.memÐf, L.memÐvf,  L.memŁvf
        @! Ðf     = Ð(f)
        @! Ð∇f    = ∇ᵢ * Ðf
        @! Ł∇f    = Ł(Ð∇f)
        @! df_dt  = p' * Ł∇f

        # dδϕ/dt
        δfᵀ_∇f, M⁻¹_δfᵀ_∇f, Ð_M⁻¹_δfᵀ_∇f = L.memŁvϕ, L.memŁvϕ, L.memÐvϕ
        @! δfᵀ_∇f       = tuple_adjoint(Łδf) * Ł∇f
        @! M⁻¹_δfᵀ_∇f   = M⁻¹ * δfᵀ_∇f
        @! Ð_M⁻¹_δfᵀ_∇f = Ð(M⁻¹_δfᵀ_∇f)
        @! dδϕ_dt       = -∇ⁱ' * Ð_M⁻¹_δfᵀ_∇f
        memÐϕ = L.memÐϕ
        for i=1:2, j=1:2
            dδϕ_dt .+= (@! memÐϕ = ∇ⁱ[i]' * (@! memÐϕ = ∇ᵢ[j]' * (@! memÐϕ = Ð(@. L.memŁϕ = t * p[j].diag * M⁻¹_δfᵀ_∇f[i]))))
        end

        FieldTuple(df_dt, dδf_dt, dδϕ_dt)

    end

    return (v!, FieldTuple(Ł(f₀), Ð(δf₀), batch_promote!(L.memÐϕ,Ð(zero(getϕ(L))))))

end

# adapting storage
adapt_structure(storage, Lϕ::LenseFlow{I,t₀,t₁}) where {I<:ODESolver,t₀,t₁} = LenseFlow{I,t₀,t₁}(adapt(storage,Lϕ.ϕ))
adapt_structure(storage, Lϕ::CachedLenseFlow{N,t₀,t₁}) where {N,t₀,t₁} = cache(LenseFlow{RK4Solver{N},t₀,t₁}(adapt(storage,Lϕ.ϕ[])), adapt(storage,Lϕ.memŁf))


text:
module Load2_module
    using Interpolations

    struct Load2_module_param
        LAMBDA::Float64
        THRESHOLD::Float64
        K::Float64
    end

    mutable struct Load2_module_variables
        x::Array{Float64}
        yitp::Union{Interpolations.GriddedInterpolation{Float64, 1, Float64, Interpolations.Gridded{Interpolations.Linear}, Tuple{Vector{Float64}}}, Nothing}
    end
end


text:
import Pkg, Downloads, ProgressMeter
using PyCall, Conda, JSON, MacroTools
import MacroTools: splitdef, combinedef

## PYTHON MODULES
const Vosk        = PyNULL()
const Sounddevice = PyNULL()
const Wave        = PyNULL()
const Zipfile     = PyNULL()
const Pynput      = PyNULL()
const Key         = PyNULL()

function __init__()
    do_restart = false
    ENV["CONDA_JL_USE_MINIFORGE"] = "1"                             # Force usage of miniforge
    if !Conda.USE_MINIFORGE
        @info "Rebuilding Conda.jl for using Miniforge..."
        Pkg.build("Conda")
        do_restart = true
    end
    ENV["PYTHON"] = ""                                              # Force PyCall to use Conda.jl
    if !any(startswith.(PyCall.python, DEPOT_PATH))                 # Rebuild of PyCall if it has not been built with Conda.jl
        @info "Rebuilding PyCall for using Julia Conda.jl..."
        Pkg.build("PyCall")
        do_restart = true
    end
    if do_restart
        @info "...rebuild completed. Restart Julia and JustSayIt."
        exit()
    end
    copy!(Vosk,        pyimport_pip("vosk"))
    copy!(Sounddevice, pyimport_pip("sounddevice"; dependency="portaudio"))
    copy!(Wave,        pyimport("wave"))
    copy!(Zipfile,     pyimport("zipfile"))
    copy!(Pynput,      pyimport_pip("pynput"))
    copy!(Key,         Pynput.keyboard.Key)
end

## CONSTANTS

const SAMPLERATE = 44100                #[Hz]
const AUDIO_READ_MAX = 512              #[bytes]
const AUDIO_ALLOC_GRANULARITY = 1024^2  #[bytes]
const AUDIO_HISTORY_MIN = 1024^2        #[bytes]
const AUDIO_ELTYPE = Int16
const AUDIO_IN_CHANNELS = 1
const COMMAND_NAME_SLEEP = "sleep"
const COMMAND_NAME_AWAKE = "awake"
const COMMAND_ABORT = "abortus"
const VARARG_END = "terminus"
const VALID_VOICEARGS_KWARGS = Dict(:model=>String, :valid_input=>AbstractArray{String}, :valid_input_auto=>Bool, :interpret_function=>Function, :use_max_speed=>Bool, :vararg_end=>String, :vararg_max=>Integer, :vararg_timeout=>AbstractFloat, :ignore_unknown=>Bool)
const DEFAULT_MODEL_NAME = "default"
const DEFAULT_RECORDER_ID = "default"
const DEFAULT_READER_ID = "default"
const TYPE_MODEL_NAME = "type"
const COMMAND_RECOGNIZER_ID = "" # NOTE: This is a safe ID as it cannot be taken by any model (raises error).
const NOISES_ENGLISH = ["huh"]
const DIGITS_ENGLISH = Dict("zero"=>"0", "one"=>"1", "two"=>"2", "three"=>"3", "four"=>"4", "five"=>"5", "six"=>"6", "seven"=>"7", "eight"=>"8", "nine"=>"9", "dot"=>".", "space"=>" ")
const ALPHABET_ENGLISH = Dict("a"=>"a", "b"=>"b", "c"=>"c", "d"=>"d", "e"=>"e", "f"=>"f", "g"=>"g", "h"=>"h", "i"=>"i", "j"=>"j", "k"=>"k", "l"=>"l", "m"=>"m", "n"=>"n", "o"=>"o", "p"=>"p", "q"=>"q", "r"=>"r", "s"=>"s", "t"=>"t", "u"=>"u", "v"=>"v", "w"=>"w", "x"=>"x", "y"=>"y", "z"=>"z", "space"=>" ")

const UNKNOWN_TOKEN = "[unk]"
const PyKey = Union{Char, PyObject}

@static if Sys.iswindows()
    const JSI_DATA        = joinpath(ENV["APPDATA"], "JustSayIt")
    const MODELDIR_PREFIX = joinpath(JSI_DATA, "models")
    const CONFIG_PREFIX   = joinpath(JSI_DATA, "config")
elseif Sys.isapple()
    const JSI_DATA        = joinpath(homedir(), "Library", "Application Support", "JustSayIt")
    const MODELDIR_PREFIX = joinpath(JSI_DATA, "models")
    const CONFIG_PREFIX   = joinpath(JSI_DATA, "config")
else
    const MODELDIR_PREFIX = joinpath(homedir(), ".local", "share", "JustSayIt", "models")
    const CONFIG_PREFIX   = joinpath(homedir(), ".config", "JustSayIt")
end

const DEFAULT_MODELDIRS = Dict(DEFAULT_MODEL_NAME => joinpath(MODELDIR_PREFIX, "vosk-model-small-en-us-0.15"),
                               TYPE_MODEL_NAME    => joinpath(MODELDIR_PREFIX, "vosk-model-en-us-daanzu-20200905"))
const DEFAULT_NOISES    = Dict(DEFAULT_MODEL_NAME => NOISES_ENGLISH,
                               TYPE_MODEL_NAME    => NOISES_ENGLISH)

DEFAULT_MODEL_REPO                 = "https://alphacephei.com/vosk/models"
DEFAULT_ENGLISH_MODEL_ARCHIVE      = "vosk-model-small-en-us-0.15.zip"
DEFAULT_ENGLISH_TYPE_MODEL_ARCHIVE = "vosk-model-en-us-daanzu-20200905.zip"

## FUNCTIONS

function pyimport_pip(modulename::AbstractString; dependency::AbstractString="", channel::AbstractString="conda-forge")
    try
        pyimport(modulename)
    catch e
        if isa(e, PyCall.PyError)
            Conda.pip_interop(true)
            Conda.pip("install", modulename)
            try
                pyimport(modulename)
            catch e
                if isa(e, PyCall.PyError) && (dependency != "") # If the module import still failed after installation, try installing the dependency with Conda first.
                    Conda.pip("uninstall --yes", modulename)
                    Conda.add(dependency; channel=channel)
                    Conda.pip("install", modulename)
                    pyimport(modulename)
                else
                    rethrow(e)
                end
            end
        else
            rethrow(e)
        end
    end
end

let
    global tic, toc
    t0::Union{Float64,Nothing} = nothing

    tic()::Float64            = ( t0 = time() )
    toc()::Float64            = ( time() - t0 )
    toc(t1::Float64)::Float64 = ( time() - t1 )
end

function pretty_dict_string(dict::Dict{String,<:Any})
    key_length_max = maximum(length.(keys(dict)))
    return join([map(sort([keys(dict)...])) do x
                    join((x, dict[x]), " "^(key_length_max+1-length(x)) * "=> ")
                end...
                ], "\n")
end

pretty_cmd_string(cmd::Union{Tuple,NTuple})   = join(map(pretty_cmd_string, cmd), " + ")
pretty_cmd_string(cmd::PyObject)              = cmd.name
pretty_cmd_string(cmd)                        = string(cmd)


text:
include("10_two_site_circuit_optimization.jl")

println("
#######################################################
# Tutorial 11                                         #
#                                                     #
# 2-site state preparation via gradient               #
# optimization                                        #
#######################################################
")

#
# State preparation
#

# Starting/reference state `|0⟩ = |Z+Z+⟩`
ψ⁰ = Zp1 * Zp2

# Target state `|ψ⟩ = (|Z+Z+⟩ + |Z-Z-⟩) / √2`
ψ = (ZpZp + ZmZm) / √2

# Find `U(θ)` that minimizes
#
# F(θ) = -|⟨ψ|U(θ)|0⟩|² = -|⟨ψ|θ⟩|²
function F(θ)
  # |θ⟩ = U(θ)|0⟩
  ψᶿ = apply(U(θ, i1, i2), ψ⁰)

  # -|⟨ψ|θ⟩|²
  return -abs(inner(ψ, ψᶿ))^2
end

println("Circuit optimization")
θ⁰ = [0, 0, 0, 0]
∂F(θ) = gradient(F, θ)[1]
θ = minimize(F, ∂F, θ⁰; nsteps=50, γ=0.1)

@show θ
@show F(θ⁰), norm(∂F(θ⁰))
@show F(θ), norm(∂F(θ))


text:
"""
    with_period(kernel::Kernel, period::Real)

Construct a transformed kernel with `period`.

# Examples

```jldoctest
julia> kernel = with_period(SqExponentialKernel(), π/2);

julia> x = rand();

julia> y = rand();

julia> kernel(x, y) ≈ (SqExponentialKernel() ∘ PeriodicTransform(2/π))(x, y)
true
```
"""
function with_period(kernel::Kernel, period::Real)
    return kernel ∘ PeriodicTransform(inv(period))
end


text:
a = [11 12 13 14; 21 22 23 24; 31 32 33 34]

@assert a[2,3] == 23                     # element in row 2, col 3
@assert a[2,2:4] == [22,23,24]           # row 2, cols 2 through 4
@assert a[1:3,3] == [13; 23; 33]         # rows 1 through 3 of col 3
@assert a[3,1:end] == [31,32,33,34]      # row 3, all elements
@assert a[3,:] == [31,32,33,34]          # row 3, all elements
@assert a[2,[1;3;4]] == [21,23,24]       # row 2, cols 1, 3, 4
@assert a[[1;3],[1;4]] == [11 14; 31 34] # very disjointed subarray
@assert a[7] == 13                       # enumerates by columns

max_stars_repo_name: idrougge/ple

max_issues_repo_name: idrougge/ple

max_forks_repo_name: idrougge/ple


text:
using ..HSA: hsa_executable_t, hsa_executable_symbol_t, BrigModuleHeader,
             hsa_isa_t, brig, finalize_brig

"Stores the necessary information about a kernel to later invoke it."
type KernelInfo
    func::Function
    argtypes::Array{DataType,1}
    executable::hsa_executable_t
    kernel_symbol::hsa_executable_symbol_t
    kernel_object
    kernarg_size
    group_size
    private_size

    function KernelInfo(executable, kernel::Function, types)
        kernel_symbol = find_kernel_symbol(executable, kernel)
        kernel_object = HSA.executable_symbol_info_kernel_object(kernel_symbol)
        kernarg_segment_size = HSA.executable_symbol_info_kernel_kernarg_segment_size(kernel_symbol)
        group_segment_size = HSA.executable_symbol_info_kernel_group_segment_size(kernel_symbol)
        private_segment_size = HSA.executable_symbol_info_kernel_private_segment_size(kernel_symbol)

        debug_print("""KernelInfo:
           $(string(kernel))($(join(types, ",")))
           kernel_object: $kernel_object
           kernarg: $kernarg_segment_size
           group: $group_segment_size
           private: $private_segment_size""")

        return new(
            kernel,
            types,
            executable,
            kernel_symbol,
            kernel_object,
            kernarg_segment_size,
            group_segment_size,
            private_segment_size)
    end
end

type KernelBlob
    brig::Ptr{BrigModuleHeader}
    executable_by_isa
    function KernelBlob(brig)
        new(brig, Dict{AbstractString, KernelInfo}())
    end
end

const kernel_cache = Dict{Tuple{Function, Array{DataType, 1}}, KernelBlob}()

function find_kernel_symbol(executable, kernel::Function)
    kernel_name = "&" * match(r"[^.]+$", string(kernel)).match

    symbols = HSA.symbols(executable)
    for (symbol, symbol_name) in symbols
        if symbol_name == kernel_name
            return symbol
        end
    end

    warn("known symbols: $(join(map(x -> x[2], symbols))) ")
    error("symbol for kernel $kernel_name not found in given executable")
end

export build_kernel
function build_kernel(agent, kernel, types)
    global kernel_cache

    debug_print("build_kernel: Kernel $(string(kernel))($(join(types,','))) for Agent '$(HSA.agent_info_name(agent))'")

    blob = get!(kernel_cache, (kernel, types)) do
        debug_print("build_kernel: build BRIG")
        KernelBlob(brig(kernel, types))
    end

    isa = HSA.agent_info_isa(agent)
    isa_name = HSA.isa_info_name(isa)
    debug_print("build_kernel: for ISA '$(isa_name)' ($(isa.handle))")

    kernel_info = get!(blob.executable_by_isa, isa_name) do
        debug_print("build_kernel: build Executable")
        code_object = finalize_brig(blob.brig, isa)

        executable = Executable()
        HSA.executable_load_code_object(executable, agent, code_object)

        HSA.executable_freeze(executable)

        return KernelInfo(executable, kernel, types)
    end

    return kernel_info
end


text:
# https://github.com/jheinen/GR.jl/issues/278#issuecomment-587090846
ENV["GKSwstype"] = "nul"

# get path of docs folders
cdir = dirname(@__FILE__);
ldir = joinpath(cdir,"literate");
tdir = joinpath(cdir,"src","literate");

# add look up paths
push!(LOAD_PATH, joinpath(cdir, ".."));

using Test, Documenter, DocumenterCitations;
using Literate;

using MagMechFEM_Matlab2Julia;

# check if run by CI
CIflag = get(ENV,"CI","") != "";

# directory filtering function
filterdir(ext,path) = filter(x->(contains(x,ext) && isfile(joinpath(path,x))),
                             readdir(path));

# convert  literate files to markdown
for f in filterdir(".jl",ldir)
  if CIflag
    Literate.markdown(joinpath(ldir,f),tdir);
  else
    Literate.markdown(joinpath(ldir,f),tdir,
                      repo_root_url="../../..");
  end
end

bib = CitationBibliography(joinpath(cdir, "bibliography.bib"))

htmlwriter = Documenter.HTML(
              collapselevel = 1,
              prettyurls = CIflag,
              assets = [asset("assets/logo.ico",class=:ico,islocal=true)]
             );

pages = ["Home"       => "index.md",
         "Packages"   => [ "NURBS Toolbox" => "NURBStoolbox.md", ],
         "API"        => [ "NURBS Toolbox" => "api_NURBStoolbox.md", ],
         "Examples"   => [ "NURBS Toolbox" => "literate/ex_NURBStoolbox.md", ],
         "References" => "references.md",
        ];

makedocs(bib,
         sitename = "MagMechFEM_Matlab2Julia",
         authors = "J. A. Duffek",
         modules = [MagMechFEM_Matlab2Julia],
         format = htmlwriter,
         doctest = false,
         clean = true,
         pages = pages
        );

if CIflag
  deploydocs(
      repo = "github.com/LazyScholar/MagMechFEM_Matlab2Julia.git",
      push_preview = true,
      forcepush = true,
      devbranch = "main",
      branch = "gh-pages"
     );
else
  println("documentation at: file:///" * joinpath(cdir,"build/index.html"));
end


text:
@testset "entities" begin
    # g, train_X, train_y = traindata(Entities(), :bgs)
    # @test typeof(g) == Dict{Any,Any}
    # @test typeof(train_X) == SparseMatrixCSC{Float32,Int64}
    # @test size(train_X) == (140, 1433)
    # @test typeof(train_y) == SparseMatrixCSC{Int32,Int64}
    # @test size(train_y) == (140, 7)

    # g, test_X, test_y = testdata(Entities(), :bgs)
    # @test typeof(g) == Dict{Any,Any}
    # @test typeof(test_X) == SparseMatrixCSC{Float32,Int64}
    # @test size(test_X) == (1000, 1433)
    # @test typeof(test_y) == SparseMatrixCSC{Int32,Int64}
    # @test size(test_y) == (1000, 7)
end


text:
using HDF5
using MultivariateStats
using Random
using CSV
using NPZ
using DPMMSubClusters
import YAML

# Define data params:
cfg = YAML.load_file("/home/vlad/projects/DPMMSubClusters.jl/benchmarks/configs/basic.yaml")

Ns = cfg["data"]["Ns"]
Ds = cfg["data"]["Ds"]
Ks = cfg["data"]["Ks"]
SEED = cfg["data"]["seed"]
DS_PATH = cfg["data"]["path"]
datavar = 10

# Define path:
mkpath(DS_PATH)

# utils function

function create_gaussian_data(path, seed, N, D, K, data_var)
	Random.seed!(seed) #generate same data
    x,y,clusters =  DPMMSubClusters.generate_gaussian_data(N, D, K, data_var)

    file_name = path * "/gaussian_D=" * string(D) *"_K=" * string(K)
    file_name_h5 = file_name*".h5"

    ispath(file_name_h5) && rm(file_name_h5)
    h5write(file_name_h5, "data", x)
    h5write(file_name_h5, "labels", y)

	npzwrite(file_name*"_x.npy", x')
	npzwrite(file_name*"_y.npy", y)

end

# generate and save data:

data_gen_path = DS_PATH * "/seed_" * string(SEED)
mkpath(data_gen_path)

# iterate over number of points:
for l in 1:length(Ns)
    N = Ns[l]

    # iterate over dimensions:
    for i in 1:length(Ds)
        D = Ds[i]

        # iterate over number of clusters
        for j in 1:length(Ks)
            K = Ks[j]

            create_gaussian_data(data_gen_path, SEED, N, D, K, datavar*j*2)

        end
    end
end


text:
import Base: size, length, iterate, getindex, setindex!, lastindex
using Base: @propagate_inbounds

import Adapt
using OffsetArrays

import Oceananigans.Architectures: architecture
import Oceananigans.Utils: datatuple
import Oceananigans.Grids: total_size, topology

using Oceananigans.Architectures
using Oceananigans.Grids
using Oceananigans.Utils

@hascuda using CuArrays

"""
    AbstractField{X, Y, Z, A, G}

Abstract supertype for fields located at `(X, Y, Z)` with data stored in a container
of type `A`. The field is defined on a grid `G`.
"""
abstract type AbstractField{X, Y, Z, A, G} end

"""
    Field{X, Y, Z, A, G, B} <: AbstractField{X, Y, Z, A, G}

A field defined at the location (`X`, `Y`, `Z`), each of which can be either `Cell`
or `Face`, and with data stored in a container of type `A` (typically an array).
The field is defined on a grid `G` and has field boundary conditions `B`.
"""
struct Field{X, Y, Z, A, G, B} <: AbstractField{X, Y, Z, A, G}
                   data :: A
                   grid :: G
    boundary_conditions :: B

    function Field{X, Y, Z}(data, grid, bcs) where {X, Y, Z}
        Tx, Ty, Tz = total_size((X, Y, Z), grid)

        if size(data) != (Tx, Ty, Tz)
            e = "Cannot construct field at ($X, $Y, $Z) with size(data)=$(size(data)). " *
                "`data` must have size ($Tx, $Ty, $Tz)."
            throw(ArgumentError(e))
        end

        return new{X, Y, Z, typeof(data), typeof(grid), typeof(bcs)}(data, grid, bcs)
    end
end

"""
    Field(L::Tuple, arch, grid, bcs, [data=zeros(arch, grid)])

Construct a `Field` on some architecture `arch` and a `grid` with some `data`.
The field's location is defined by a tuple `L` of length 3 whose elements are
`Cell` or `Face` and has field boundary conditions `bcs`.
"""
function Field(L::Tuple, arch, grid, bcs,
               data=zeros(eltype(grid), arch, grid, (typeof(L[1]), typeof(L[2]), typeof(L[3]))))

    return Field{typeof(L[1]), typeof(L[2]), typeof(L[3])}(data, grid, bcs)
end

function Field(L::NTuple{3, DataType}, arch, grid, bcs,
               data=zeros(eltype(grid), arch, grid, (L[1], L[2], L[3])))

    return Field{L[1], L[2], L[3]}(data, grid, bcs)
end

"""
    Field(X, Y, Z, arch, grid, [data=zeros(arch, grid)], bcs)

Construct a `Field` on some architecture `arch` and a `grid` with some `data`.
The field's location is defined by `X`, `Y`, `Z` where each is either `Cell` or `Face`
and has field boundary conditions `bcs`.
"""
Field(X, Y, Z, arch, grid, bcs, data=zeros(eltype(grid), arch, grid, (X, Y, Z))) =
    Field((X, Y, Z), arch, grid, bcs, data)

"""
    CellField([FT=eltype(grid)], arch::AbstractArchitecture, grid, bcs=TracerBoundaryConditions(grid),
              data=zeros(FT, arch, grid, (Cell, Cell, Cell)))

Return a `Field{Cell, Cell, Cell}` on architecture `arch` and `grid` containing `data`
with field boundary conditions `bcs`.
"""
function CellField(arch::AbstractArchitecture, grid,
                   bcs=TracerBoundaryConditions(grid),
                   data=zeros(eltype(grid), arch, grid, (Cell, Cell, Cell)))

    return Field(Cell, Cell, Cell, arch, grid, bcs, data)
end

"""
    XFaceField([FT=eltype(grid)], arch::AbstractArchitecture, grid, bcs=UVelocityBoundaryConditions(grid),
              data=zeros(FT, arch, grid, (Face, Cell, Cell)))

Return a `Field{Face, Cell, Cell}` on architecture `arch` and `grid` containing `data`
with field boundary conditions `bcs`.
"""
function XFaceField(arch::AbstractArchitecture, grid,
                    bcs=UVelocityBoundaryConditions(grid),
                    data=zeros(eltype(grid), arch, grid, (Face, Cell, Cell)))

    return Field(Face, Cell, Cell, arch, grid, bcs, data)
end

"""
    YFaceField([FT=eltype(grid)], arch::AbstractArchitecture, grid, bcs=VVelocityBoundaryConditions(grid),
              data=zeros(FT, arch, grid, (Cell, Face, Cell)))

Return a `Field{Cell, Face, Cell}` on architecture `arch` and `grid` containing `data`
with field boundary conditions `bcs`.
"""
function YFaceField(arch::AbstractArchitecture, grid,
                    bcs=VVelocityBoundaryConditions(grid),
                    data=zeros(eltype(grid), arch, grid, (Cell, Face, Cell)))

    return Field(Cell, Face, Cell, arch, grid, bcs, data)
end

"""
    ZFaceField([FT=eltype(grid)], arch::AbstractArchitecture, grid, bcs=WVelocityBoundaryConditions(grid),
              data=zeros(FT, arch, grid, (Cell, Cell, Face))

Return a `Field{Cell, Cell, Face}` on architecture `arch` and `grid` containing `data`
with field boundary conditions `bcs`.
"""
function ZFaceField(arch::AbstractArchitecture, grid,
                    bcs=WVelocityBoundaryConditions(grid),
                    data=zeros(eltype(grid), arch, grid, (Cell, Cell, Face)))

    return Field(Cell, Cell, Face, arch, grid, bcs, data)
end

CellField(FT::DataType, arch, grid, bcs=TracerBoundaryConditions(grid)) =
    CellField(arch, grid, bcs, zeros(FT, arch, grid, (Cell, Cell, Cell)))

XFaceField(FT::DataType, arch, grid, bcs=UVelocityBoundaryConditions(grid)) =
    XFaceField(arch, grid, bcs, zeros(FT, arch, grid, (Face, Cell, Cell)))

YFaceField(FT::DataType, arch, grid, bcs=VVelocityBoundaryConditions(grid)) =
    YFaceField(arch, grid, bcs, zeros(FT, arch, grid, (Cell, Face, Cell)))

ZFaceField(FT::DataType, arch, grid, bcs=WVelocityBoundaryConditions(grid)) =
    ZFaceField(arch, grid, bcs, zeros(FT, arch, grid, (Cell, Cell, Face)))

#####
##### Functions for querying fields
#####

location(a) = nothing

"Returns the location `(X, Y, Z)` of an `AbstractField{X, Y, Z}`."
location(::AbstractField{X, Y, Z}) where {X, Y, Z} = (X, Y, Z) # note no instantiation

x_location(::AbstractField{X}) where X               = X
y_location(::AbstractField{X, Y}) where {X, Y}       = Y
z_location(::AbstractField{X, Y, Z}) where {X, Y, Z} = Z

"Returns the architecture where the field data `f.data` is stored."
architecture(f::Field) = architecture(f.data)
architecture(o::OffsetArray) = architecture(o.parent)

"Returns the length of a field's `data`."
@inline length(f::Field) = length(f.data)

"Returns the topology of a fields' `grid`."
@inline topology(f, args...) = topology(f.grid, args...)

"""
Returns the length of a field located at `Cell` centers along a grid
dimension of length `N` and with halo points `H`.
"""
dimension_length(loc, topo, N, H=0) = N + 2H

"""
Returns the length of a field located at cell `Face`s along a grid
dimension of length `N` and with halo points `H`.
"""
dimension_length(::Type{Face}, ::Bounded, N, H=0) = N + 1 + 2H

"""
    size(loc, grid)

Returns the size of a field at `loc` on `grid`.
This is a 3-tuple of integers corresponding to the number of interior nodes
of `f` along `x, y, z`.
"""
@inline size(loc, grid) = (dimension_length(loc[1], topology(grid, 1), grid.Nx),
                           dimension_length(loc[2], topology(grid, 2), grid.Ny),
                           dimension_length(loc[3], topology(grid, 3), grid.Nz))

"""
    size(f::AbstractField{X, Y, Z}) where {X, Y, Z}

Returns the size of an `AbstractField{X, Y, Z}` located at `X, Y, Z`.
This is a 3-tuple of integers corresponding to the number of interior nodes
of `f` along `x, y, z`.
"""
@inline size(f::AbstractField) = size(location(f), f.grid)

"""
    total_size(loc, grid)

Returns the "total" size of a field at `loc` on `grid`.
This is a 3-tuple of integers corresponding to the number of grid points
contained by `f` along `x, y, z`.
"""
@inline total_size(loc, grid) = (dimension_length(loc[1], topology(grid, 1), grid.Nx, grid.Hx),
                                 dimension_length(loc[2], topology(grid, 2), grid.Ny, grid.Hy),
                                 dimension_length(loc[3], topology(grid, 3), grid.Nz, grid.Hz))

@inline total_size(f::AbstractField) = total_size(location(f), f.grid)

@propagate_inbounds getindex(f::Field, inds...) = @inbounds getindex(f.data, inds...)
@propagate_inbounds setindex!(f::Field, v, inds...) = @inbounds setindex!(f.data, v, inds...)

@inline lastindex(f::Field) = lastindex(f.data)
@inline lastindex(f::Field, dim) = lastindex(f.data, dim)

"Returns `f.data` for `f::Field` or `f` for `f::AbstractArray."
@inline data(a) = a # fallback
@inline data(f::Field) = f.data

# Endpoint for recursive `datatuple` function:
@inline datatuple(obj::AbstractField) = data(obj)

"Returns `f.data.parent` for `f::Field`."
@inline Base.parent(f::Field) = f.data.parent

@inline interior_indices(loc, topo, N) = 1:N
@inline interior_indices(::Type{Face}, ::Bounded, N) = 1:N+1

@inline interior_parent_indices(loc, topo, N, H) = 1+H:N+H
@inline interior_parent_indices(::Type{Face}, ::Bounded, N, H) = 1+H:N+1+H

"Returns a view of `f` that excludes halo points."
@inline interior(f::Field{X, Y, Z}) where {X, Y, Z} = view(f.data, interior_indices(X, topology(f, 1), f.grid.Nx),
                                                                   interior_indices(Y, topology(f, 2), f.grid.Ny),
                                                                   interior_indices(Z, topology(f, 3), f.grid.Nz))

"Returns a reference (not a view) to the interior points of `field.data.parent.`"
@inline interiorparent(f::Field{X, Y, Z}) where {X, Y, Z} =
    @inbounds f.data.parent[interior_parent_indices(X, topology(f, 1), f.grid.Nx, f.grid.Hx),
                            interior_parent_indices(Y, topology(f, 2), f.grid.Ny, f.grid.Hy),
                            interior_parent_indices(Z, topology(f, 3), f.grid.Nz, f.grid.Hz)]

iterate(f::Field, state=1) = iterate(f.data, state)

@inline xnode(::Type{Cell}, i, grid) = @inbounds grid.xC[i]
@inline xnode(::Type{Face}, i, grid) = @inbounds grid.xF[i]

@inline ynode(::Type{Cell}, j, grid) = @inbounds grid.yC[j]
@inline ynode(::Type{Face}, j, grid) = @inbounds grid.yF[j]

@inline znode(::Type{Cell}, k, grid) = @inbounds grid.zC[k]
@inline znode(::Type{Face}, k, grid) = @inbounds grid.zF[k]

@inline xnode(i, ψ::Field{X, Y, Z}) where {X, Y, Z} = xnode(X, i, ψ.grid)
@inline ynode(j, ψ::Field{X, Y, Z}) where {X, Y, Z} = ynode(Y, j, ψ.grid)
@inline znode(k, ψ::Field{X, Y, Z}) where {X, Y, Z} = znode(Z, k, ψ.grid)

# Dispatch insanity
xnodes(::Type{Cell}, topo, grid) = reshape(grid.xC, grid.Nx, 1, 1)
ynodes(::Type{Cell}, topo, grid) = reshape(grid.yC, 1, grid.Ny, 1)
znodes(::Type{Cell}, topo, grid) = reshape(grid.zC, 1, 1, grid.Nz)

xnodes(::Type{Face}, topo, grid) = reshape(grid.xF[1:end-1], grid.Nx, 1, 1)
ynodes(::Type{Face}, topo, grid) = reshape(grid.yF[1:end-1], 1, grid.Ny, 1)
znodes(::Type{Face}, topo, grid) = reshape(grid.zF[1:end-1], 1, 1, grid.Nz)

xnodes(::Type{Face}, ::Bounded, grid) = reshape(grid.xF, grid.Nx+1, 1, 1)
ynodes(::Type{Face}, ::Bounded, grid) = reshape(grid.yF, 1, grid.Ny+1, 1)
znodes(::Type{Face}, ::Bounded, grid) = reshape(grid.zF, 1, 1, grid.Nz+1)

xnodes(ψ::AbstractField) = xnodes(x_location(ψ), topology(ψ, 1), ψ.grid)
ynodes(ψ::AbstractField) = ynodes(y_location(ψ), topology(ψ, 2), ψ.grid)
znodes(ψ::AbstractField) = znodes(z_location(ψ), topology(ψ, 3), ψ.grid)

nodes(ψ::AbstractField) = (xnodes(ψ), ynodes(ψ), znodes(ψ))

#####
##### Creating offset arrays for field data by dispatching on architecture.
#####

"""
Return a range of indices for a field located at `Cell` centers
`along a grid dimension of length `N` and with halo points `H`.
"""
offset_indices(loc, topo, N, H=0) = 1 - H : N + H

"""
Return a range of indices for a field located at cell `Face`s
`along a grid dimension of length `N` and with halo points `H`.
"""
offset_indices(::Type{Face}, ::Bounded, N, H=0) = 1 - H : N + H + 1

"""
    OffsetArray(underlying_data, grid, loc)

Returns an `OffsetArray` that maps to `underlying_data` in memory,
with offset indices appropriate for the `data` of a field on
a `grid` of `size(grid)` and located at `loc`.
"""
function OffsetArray(underlying_data, grid, loc)
    ii = offset_indices(loc[1], topology(grid, 1), grid.Nx, grid.Hx)
    jj = offset_indices(loc[2], topology(grid, 2), grid.Ny, grid.Hy)
    kk = offset_indices(loc[3], topology(grid, 3), grid.Nz, grid.Hz)

    return OffsetArray(underlying_data, ii, jj, kk)
end

"""
    zeros([FT=Float64], ::CPU, grid, loc)

Returns an `OffsetArray` of zeros of float type `FT`, with
parent data in CPU memory and indices corresponding to a field on a
`grid` of `size(grid)` and located at `loc`.
"""
function Base.zeros(FT, ::CPU, grid, loc=(Cell, Cell, Cell))
    underlying_data = zeros(FT, dimension_length(loc[1], topology(grid, 1), grid.Nx, grid.Hx),
                                dimension_length(loc[2], topology(grid, 2), grid.Ny, grid.Hy),
                                dimension_length(loc[3], topology(grid, 3), grid.Nz, grid.Hz))

    return OffsetArray(underlying_data, grid, loc)
end

"""
    zeros([FT=Float64], ::GPU, grid, loc)

Returns an `OffsetArray` of zeros of float type `FT`, with
parent data in GPU memory and indices corresponding to a field on a `grid`
of `size(grid)` and located at `loc`.
"""
function Base.zeros(FT, ::GPU, grid, loc=(Cell, Cell, Cell))
    underlying_data = CuArray{FT}(undef, dimension_length(loc[1], topology(grid, 1), grid.Nx, grid.Hx),
                                         dimension_length(loc[2], topology(grid, 2), grid.Ny, grid.Hy),
                                         dimension_length(loc[3], topology(grid, 3), grid.Nz, grid.Hz))

    underlying_data .= 0 # Ensure data is initially 0.

    return OffsetArray(underlying_data, grid, loc)
end

# Default to type of Grid
Base.zeros(arch, grid, loc=(Cell, Cell, Cell)) = zeros(eltype(grid), arch, grid, loc)

Base.zeros(FT, ::CPU, grid, Nx, Ny, Nz) = zeros(FT, Nx, Ny, Nz)
Base.zeros(FT, ::GPU, grid, Nx, Ny, Nz) = zeros(FT, Nx, Ny, Nz) |> CuArray

Base.zeros(arch, grid, args...) = zeros(eltype(grid), args...)


text:
function decompose_cell(cell_connectivity, cell_type::Int; target::String="")
    if target == ""
        if in(cell_type, POINT_CELLS) || in(cell_type, LINE_CELLS)
            throw("Cell is already represented as points.")
        elseif in(cell_type, FACE_CELLS)
            target = "Lines"
        elseif in(cell_type, VOLUME_CELLS)
            target = "Faces"
        end
    end

    args = (cell_connectivity, target)
    if cell_type == 5
        return decompose_triangle(args...)
    elseif cell_type == 6
        return decompose_triangle_strip(args...)
    elseif cell_type == 8
        return decompose_pixel(args...)
    elseif cell_type == 9
        return decompose_quad(args...)
    elseif cell_type == 10
        return decompose_tetra(args...)
    elseif cell_type == 11
        return decompose_voxel(args...)
    elseif cell_type == 12
        return decompose_hexa(args...)
    elseif cell_type == 13
        return decompose_wedge(args...)
    elseif cell_type == 14
        return decompose_pyramid(args...)
    elseif cell_type == 22
        return decompose_quadratic_triangle(args...)
    elseif cell_type == 23
        return decompose_quadratic_quad(args...)
    elseif cell_type == 24
        return decompose_quadratic_tetra(args...)
    elseif cell_type == 25
        return decompose_quadratic_hexa(args...)
    else
        throw("Cell type $cell_type is not supported.")
    end
end

function decompose_triangle(cell_connectivity, target::String)
    if target == "Faces"
        return Vector{Int}[[cell_connectivity...]], Int[5]
    else
        return Vector{Int}[ [ cell_connectivity[1], cell_connectivity[2] ],
            [ cell_connectivity[2], cell_connectivity[3] ],
            [ cell_connectivity[3], cell_connectivity[1] ] ],
            fill(3,3)
    end
end

function decompose_triangle_strip(cell_connectivity, target::String)
    if target == "Faces"
        return Vector{Int}[ [ cell_connectivity[i], cell_connectivity[i+1], cell_connectivity[i+2] ]
            for i in 1:length(cell_connectivity)-2], fill(5, length(cell_connectivity)-2)
    else
        return append!(
            Vector{Int}[ [ cell_connectivity[i], cell_connectivity[i+1] ]
                for i in 1:length(cell_connectivity)-1 ],
            Vector{Int}[ [ cell_connectivity[i], cell_connectivity[i+2] ]
                for i in 1:length(cell_connectivity)-2]),
            fill(3, 2*length(cell_connectivity)-3)
    end
end

function decompose_pixel(cell_connectivity, target::String)
    if target == "Faces"
        return [cell_connectivity], [8]
    else
        return Vector{Int}[ [ cell_connectivity[1], cell_connectivity[2] ],
            [ cell_connectivity[2], cell_connectivity[4] ],
            [ cell_connectivity[4], cell_connectivity[3] ],
            [ cell_connectivity[3], cell_connectivity[1] ] ],
            fill(3, 4)
    end
end

function decompose_quad(cell_connectivity, target::String)
    if target == "Faces"
        return [cell_connectivity], [9]
    else
        return Vector{Int}[ [ cell_connectivity[1], cell_connectivity[2] ],
            [ cell_connectivity[2], cell_connectivity[3] ],
            [ cell_connectivity[3], cell_connectivity[4] ],
            [ cell_connectivity[4], cell_connectivity[1] ] ],
            fill(3, 4)
    end
end

function decompose_tetra(cell_connectivity, target::String)
    if target == "Faces"
        return Vector{Int}[ [ cell_connectivity[1], cell_connectivity[2], cell_connectivity[3] ],
            [ cell_connectivity[1], cell_connectivity[2], cell_connectivity[4] ],
            [ cell_connectivity[2], cell_connectivity[3], cell_connectivity[4] ],
            [ cell_connectivity[2], cell_connectivity[1], cell_connectivity[4] ] ],
            fill(5, 4)
    else
        return Vector{Int}[ [ cell_connectivity[1], cell_connectivity[2] ],
            [ cell_connectivity[2], cell_connectivity[3] ],
            [ cell_connectivity[3], cell_connectivity[1] ],
            [ cell_connectivity[1], cell_connectivity[4] ],
            [ cell_connectivity[2], cell_connectivity[4] ],
            [ cell_connectivity[3], cell_connectivity[4] ] ],
            fill(3, 6)
    end
end

function decompose_voxel(cell_connectivity, target::String)
    if target == "Faces"
        return Vector{Int}[ [ cell_connectivity[1], cell_connectivity[2], cell_connectivity[4], cell_connectivity[3] ],
            [ cell_connectivity[5], cell_connectivity[6], cell_connectivity[8], cell_connectivity[7] ],
            [ cell_connectivity[1], cell_connectivity[2], cell_connectivity[6], cell_connectivity[5] ],
            [ cell_connectivity[2], cell_connectivity[4], cell_connectivity[8], cell_connectivity[6] ],
            [ cell_connectivity[4], cell_connectivity[3], cell_connectivity[7], cell_connectivity[8] ],
            [ cell_connectivity[3], cell_connectivity[1], cell_connectivity[5], cell_connectivity[7] ] ],
            fill(9, 6)
    else
        return Vector{Int}[ [ cell_connectivity[1], cell_connectivity[2] ],
            [ cell_connectivity[2], cell_connectivity[4] ],
            [ cell_connectivity[4], cell_connectivity[3] ],
            [ cell_connectivity[3], cell_connectivity[1] ],
            [ cell_connectivity[5], cell_connectivity[6] ],
            [ cell_connectivity[6], cell_connectivity[8] ],
            [ cell_connectivity[8], cell_connectivity[7] ],
            [ cell_connectivity[7], cell_connectivity[5] ],
            [ cell_connectivity[1], cell_connectivity[5] ],
            [ cell_connectivity[2], cell_connectivity[6] ],
            [ cell_connectivity[4], cell_connectivity[8] ],
            [ cell_connectivity[3], cell_connectivity[7] ] ],
            fill(3, 12)
    end
end

function decompose_hexa(cell_connectivity, target::String)
    if target == "Faces"
        return Vector{Int}[ [ cell_connectivity[1], cell_connectivity[2], cell_connectivity[3], cell_connectivity[4] ],
            [ cell_connectivity[5], cell_connectivity[6], cell_connectivity[7], cell_connectivity[8] ],
            [ cell_connectivity[1], cell_connectivity[2], cell_connectivity[6], cell_connectivity[5] ],
            [ cell_connectivity[2], cell_connectivity[3], cell_connectivity[7], cell_connectivity[6] ],
            [ cell_connectivity[3], cell_connectivity[4], cell_connectivity[8], cell_connectivity[7] ],
            [ cell_connectivity[4], cell_connectivity[1], cell_connectivity[5], cell_connectivity[8] ] ],
            fill(9, 6)
    else
        return Vector{Int}[ [ cell_connectivity[1], cell_connectivity[2] ],
            [ cell_connectivity[2], cell_connectivity[3] ],
            [ cell_connectivity[3], cell_connectivity[4] ],
            [ cell_connectivity[4], cell_connectivity[1] ],
            [ cell_connectivity[5], cell_connectivity[6] ],
            [ cell_connectivity[6], cell_connectivity[7] ],
            [ cell_connectivity[7], cell_connectivity[8] ],
            [ cell_connectivity[8], cell_connectivity[5] ],
            [ cell_connectivity[1], cell_connectivity[5] ],
            [ cell_connectivity[2], cell_connectivity[6] ],
            [ cell_connectivity[3], cell_connectivity[7] ],
            [ cell_connectivity[4], cell_connectivity[8] ] ],
            fill(3, 12)
    end
end

function decompose_wedge(cell_connectivity, target::String)
    if target == "Faces"
        return Vector{Int}[ [ cell_connectivity[1], cell_connectivity[2], cell_connectivity[5], cell_connectivity[4] ],
            [ cell_connectivity[1], cell_connectivity[2], cell_connectivity[3] ],
            [ cell_connectivity[2], cell_connectivity[5], cell_connectivity[6], cell_connectivity[3] ],
            [ cell_connectivity[5], cell_connectivity[4], cell_connectivity[6] ],
            [ cell_connectivity[4], cell_connectivity[1], cell_connectivity[3], cell_connectivity[6] ] ],
            [9, 5, 9, 5, 9]
    else
        return Vector{Int}[ [ cell_connectivity[1], cell_connectivity[2] ],
            [ cell_connectivity[2], cell_connectivity[5] ],
            [ cell_connectivity[5], cell_connectivity[4] ],
            [ cell_connectivity[4], cell_connectivity[1] ],
            [ cell_connectivity[1], cell_connectivity[3] ],
            [ cell_connectivity[2], cell_connectivity[3] ],
            [ cell_connectivity[5], cell_connectivity[6] ],
            [ cell_connectivity[4], cell_connectivity[6] ],
            [ cell_connectivity[3], cell_connectivity[6] ] ],
            fill(3, 9)
    end
end

function decompose_pyramid(cell_connectivity, target::String)
    if target == "Faces"
        return Vector{Int}[ [ cell_connectivity[1], cell_connectivity[2], cell_connectivity[3], cell_connectivity[4] ],
            [ cell_connectivity[1], cell_connectivity[2], cell_connectivity[5] ],
            [ cell_connectivity[2], cell_connectivity[3], cell_connectivity[5] ],
            [ cell_connectivity[3], cell_connectivity[4], cell_connectivity[5] ],
            [ cell_connectivity[4], cell_connectivity[1], cell_connectivity[5] ] ],
            Int[9, 5, 5, 5, 5]
    else
        return Vector{Int}[ [ cell_connectivity[1], cell_connectivity[2] ],
            [ cell_connectivity[2], cell_connectivity[3] ],
            [ cell_connectivity[3], cell_connectivity[4] ],
            [ cell_connectivity[4], cell_connectivity[1] ],
            [ cell_connectivity[1], cell_connectivity[5] ],
            [ cell_connectivity[2], cell_connectivity[5] ],
            [ cell_connectivity[3], cell_connectivity[5] ],
            [ cell_connectivity[4], cell_connectivity[5] ] ],
            fill(3, 8)
    end
end

function decompose_quadratic_triangle(cell_connectivity, target::String)
    if target == "Faces"
        return Vector{Int}[cell_connectivity], Int[22]
    else
        return Vector{Int}[ [ cell_connectivity[1], cell_connectivity[2], cell_connectivity[4]],
            [ cell_connectivity[2], cell_connectivity[3], cell_connectivity[5] ],
            [ cell_connectivity[3], cell_connectivity[1], cell_connectivity[6] ] ],
            fill(21, 3)
    end
end

function decompose_quadratic_quad(cell_connectivity, target::String)
    if target == "Faces"
        return Vector{Int}[cell_connectivity], Int[23]
    else
        return Vector{Int}[ [ cell_connectivity[1], cell_connectivity[2], cell_connectivity[5] ],
            [ cell_connectivity[2], cell_connectivity[3], cell_connectivity[6] ],
            [ cell_connectivity[3], cell_connectivity[4], cell_connectivity[7] ],
            [ cell_connectivity[4], cell_connectivity[1], cell_connectivity[8] ] ],
            fill(21, 4)
    end
end

function decompose_quadratic_tetra(cell_connectivity, target::String)
    if target == "Faces"
        return Vector{Int}[ [ cell_connectivity[1], cell_connectivity[2], cell_connectivity[3], cell_connectivity[5], cell_connectivity[6], cell_connectivity[7] ],
            [ cell_connectivity[1], cell_connectivity[2], cell_connectivity[4], cell_connectivity[5], cell_connectivity[9], cell_connectivity[8] ],
            [ cell_connectivity[2], cell_connectivity[3], cell_connectivity[4], cell_connectivity[6], cell_connectivity[10], cell_connectivity[9] ],
            [ cell_connectivity[3], cell_connectivity[1], cell_connectivity[4], cell_connectivity[7], cell_connectivity[8], cell_connectivity[10] ] ],
            fill(22, 4)
    else
        return Vector{Int}[ [ cell_connectivity[1], cell_connectivity[2], cell_connectivity[5] ],
            [ cell_connectivity[2], cell_connectivity[3], cell_connectivity[5] ],
            [ cell_connectivity[3], cell_connectivity[1], cell_connectivity[7] ],
            [ cell_connectivity[1], cell_connectivity[4], cell_connectivity[8] ],
            [ cell_connectivity[2], cell_connectivity[4], cell_connectivity[8] ],
            [ cell_connectivity[3], cell_connectivity[4], cell_connectivity[10] ] ],
            fill(21, 6)
    end
end

function decompose_quadratic_hexa(cell_connectivity, target::String)
    if target == "Faces"
        return Vector{Int}[ [ cell_connectivity[1], cell_connectivity[2], cell_connectivity[3], cell_connectivity[4], cell_connectivity[9], cell_connectivity[10], cell_connectivity[11], cell_connectivity[12] ],
            [ cell_connectivity[5], cell_connectivity[6], cell_connectivity[7], cell_connectivity[8], cell_connectivity[13], cell_connectivity[14], cell_connectivity[15], cell_connectivity[16] ],
            [ cell_connectivity[1], cell_connectivity[2], cell_connectivity[6], cell_connectivity[5], cell_connectivity[9], cell_connectivity[18], cell_connectivity[13], cell_connectivity[17] ],
            [ cell_connectivity[2], cell_connectivity[3], cell_connectivity[7], cell_connectivity[6], cell_connectivity[10], cell_connectivity[19], cell_connectivity[14], cell_connectivity[18] ],
            [ cell_connectivity[3], cell_connectivity[4], cell_connectivity[8], cell_connectivity[7], cell_connectivity[11], cell_connectivity[20], cell_connectivity[15], cell_connectivity[19] ],
            [ cell_connectivity[4], cell_connectivity[1], cell_connectivity[5], cell_connectivity[8], cell_connectivity[12], cell_connectivity[17], cell_connectivity[16], cell_connectivity[20] ] ],
            fill(23, 6)
    else
        return Vector{Int}[ [ cell_connectivity[1], cell_connectivity[2], cell_connectivity[9] ],
            [ cell_connectivity[2], cell_connectivity[3], cell_connectivity[10] ],
            [ cell_connectivity[3], cell_connectivity[4], cell_connectivity[11] ],
            [ cell_connectivity[4], cell_connectivity[1], cell_connectivity[12] ],
            [ cell_connectivity[5], cell_connectivity[6], cell_connectivity[13] ],
            [ cell_connectivity[6], cell_connectivity[7], cell_connectivity[14] ],
            [ cell_connectivity[7], cell_connectivity[8], cell_connectivity[15] ],
            [ cell_connectivity[8], cell_connectivity[5], cell_connectivity[16] ],
            [ cell_connectivity[1], cell_connectivity[5], cell_connectivity[17] ],
            [ cell_connectivity[2], cell_connectivity[6], cell_connectivity[18] ],
            [ cell_connectivity[3], cell_connectivity[7], cell_connectivity[19] ],
            [ cell_connectivity[4], cell_connectivity[8], cell_connectivity[20] ] ],
            fill(21, 12)
    end
end

function decompose_cell(dataset::AbstractStaticVTKData, cell_ind::Union{Int, Tuple{Int, Vararg{Int}}}; target::String="")
    decompose_cell(cell_connectivity(dataset, cell_ind), cell_type(dataset, cell_ind))
end


text:
__precompile__()

module ObsDataStructures
import Base: getindex, +
export label_treatments,
       Counterfactuals, Outcome, Covariates, ObsData,
       getindex, Γ, +

function label_treatments(W::Vector)
    treatment_levels = sort!([w for w in Set(W)])
    if length(treatment_levels) != 2
        error("Only binary treatments")
    end
    return Dict(label=>t for (label, t) in zip([true, false], treatment_levels))
end

type Counterfactuals{T<:Real}
    treated::Dict{Bool,Vector{T}}
    observed::Dict{Bool,Vector{T}}
    W::Vector{Bool}
end

function Counterfactuals{T<:Real}(treated::Dict{Bool,Vector{T}}, W::Vector{Bool})
    observed = Dict(c=>Vector{T}(length(W)) for c in [true, false])
    observed[true][W] = treated[true][W]
    observed[true][!W] = treated[false][!W]
    observed[false][W] = treated[false][W]
    observed[false][!W] = treated[true][!W]
    return Counterfactuals(treated, observed, W)
end

function getindex(F::Counterfactuals, index)
    return Counterfactuals(Dict(true=>F.treated[true][index],
                                false=>F.treated[false][index]),
                           F.W[index])
end

function Γ(F::Counterfactuals; o::Union{Bool,Void}=nothing, t::Union{Bool,Void}=nothing)
    if (o==nothing) & (t==nothing)
        error("no conditions provided")
    elseif o==nothing
        return F.treated[t]
    elseif t==nothing
        return F.observed[o]
    else
        return t ? F.observed[o][F.W] : F.observed[o][!F.W]
    end
end

function +(F::Counterfactuals, f::Counterfactuals)
    treated = Dict(t=>F.treated[t]+f.treated[t] for t in (true, false))
    observed = Dict(t=>F.observed[t]+f.observed[t] for t in (true, false))
    return Counterfactuals(treated, observed, F.W)
end

immutable Outcome{T<:Real}
    Y::Vector{T}
    W::Vector{Bool}
end

function getindex(Y::Outcome, index)
    return Outcome(Y.Y[index], Y.W[index])
end

function Γ(Y::Outcome; t::Bool=true)
    t ? Y.Y[Y.W] : Y.Y[!Y.W]
end

function -(Y::Outcome, F::Counterfactuals)
    return Outcome(Y.Y - F.observed[true], Y.W)
end

immutable Covariates{T<:Real}
    X::Matrix{T}
    W::Vector{Bool}
end

function getindex(X::Covariates, index)
    return Covariates(X.X[index,:], X.W[index])
end

function Γ(X::Covariates; t::Bool=true)
    t ? X.X[X.W,:] : X.X[!X.W,:]
end

struct ObsData
    X::Covariates
    W::Vector{Bool}
    Y::Outcome
    trts::Dict{Bool,Any}
    N::Int
    N_treated::Dict{Bool,Int}
end

function ObsData(X::Matrix, W::Vector, Y::Vector; trts=nothing)
    trts == nothing ? trts = label_treatments(W) : nothing
    W_bool = Vector{Bool}(length(W))
    for (t,w) in trts
        W_bool[W.==w] = t
    end
    return ObsData(Covariates(X,W_bool), W_bool, Outcome(Y, W_bool), trts, length(W), Dict(true=>sum(W_bool), false=>sum(!W_bool)))
end

function getindex(data::ObsData, index)
    return ObsData(data.X[index], data.W[index], data.Y[index], data.trts,
                   length(data.W), Dict(true=>sum(data.W), false=>sum(!data.W)))
end

end #module

max_stars_repo_name: naskoD/bachelorThesis

max_issues_repo_name: naskoD/bachelorThesis

max_forks_repo_name: naskoD/bachelorThesis


text:
#                _       _ _
#  ___  ___ _ __(_) __ _| (_)_______
# / __|/ _ \ '__| |/ _` | | |_  / _ \
# \__ \  __/ |  | | (_| | | |/ /  __/
# |___/\___|_|  |_|\__,_|_|_/___\___|

const null_attr = Dict{AbstractString, Any}("NULL" => true)

# general objects
function attribute_value(x)
    r = Dict{String, Any}()
    for e=fieldnames(x)
        r[string(e)] = null_or_val(getfield(x,e))
    end
    Dict{String, Any}("M" => r)
end

# TODO -- binary?
attribute_value(x :: Bool) = Dict{AbstractString, Any}("BOOL" => x)
attribute_value(x :: Real) = Dict{AbstractString, Any}("N" => string(x))
attribute_value(x :: AbstractString) = Dict{AbstractString, Any}("S" => x)

attribute_value(x :: Array) =
    Dict{AbstractString, Any}("L" => [attribute_value(e) for e=x])
attribute_value(x :: Set{T}) where {T <: Real} =
    Dict{AbstractString, Any}("NS" => [string(e) for e=x])
attribute_value(x :: Set{T}) where {T <: AbstractString} =
    Dict{AbstractString, Any}("SS" => [e for e=x])
# TODO -- n-dimensional arrays

function attribute_value(x :: Dict)
    dict = Dict{String, Any}()
    for (k,v)=x
        dict[string(k)] = null_or_val(v)
    end
    Dict{String, Any}("M" => dict)
end

null_or_val(x) = x == nothing ? null_attr : attribute_value(x)

#      _                     _       _ _
#   __| | ___  ___  ___ _ __(_) __ _| (_)_______
#  / _` |/ _ \/ __|/ _ \ '__| |/ _` | | |_  / _ \
# | (_| |  __/\__ \  __/ |  | | (_| | | |/ /  __/
#  \__,_|\___||___/\___|_|  |_|\__,_|_|_/___\___|

bool_val(val :: Bool) = val
function bool_val(val :: AbstractString)
    if lowercase(val) == "false"
        return false
    elseif lowercase(val) == "true"
        return true
    else
        error("Received non-boolean value for boolean typed data: $val")
    end
end

real_val(val :: Real) = val
real_val(val :: AbstractString) = parse(Float64, val)

function value_from_attributes(hash :: Dict)

    ks = keys(hash)
    if length(ks) != 1
        error("multiple type keys provided in DynamoDB typed JSON value")
    end

    ty = first(ks)
    val = hash[ty]

    if ty == "NULL"     ; return nothing
    elseif ty == "BOOL" ; return bool_val(val)
    elseif ty == "N"    ; return real_val(val)
    elseif ty == "S"    ; return val
    elseif ty == "L"    ; return [value_from_attributes(e) for e in val]
    elseif ty == "NS"   ; return Set{Real}([float(e) for e in val])
    elseif ty == "SS"   ; return Set{AbstractString}(val)
    elseif ty == "M"
        m = Dict{AbstractString, Any}()
        for (k, v)=val
            m[string(k)] = value_from_attributes(v)
        end
        return m
    else
        error("Unknown datatype value in DynamoDB typed JSON value: $ty")
    end
end

function value_from_attributes(ty :: Type, typed_json :: Dict)

    vals = value_from_attributes(Dict("M" => typed_json))

    if ty <: Dict
      return vals
    else
      init_vals = [vals[string(e)] for e=fieldnames(ty)]
      return ty(init_vals...)
    end

end


text:
using sb, LinearAlgebra, DelimitedFiles, StaticArrays

struct storages
    ts :: Vector{Float64}
    sz :: Vector{Float64}
    function storages()
        new(Vector{Float64}(),Vector{Float64}())
    end
end

function store!(storage, t, ops, sbm)
    append!(storage.ts, t)
    append!(storage.sz, real(tr(sb.σz*ops.ρ)))
end

const NSTEPS = 2500
const DT = 0.010
const TOL = 1.0e-11

function my_run(sbm)
    storage = storages()
    ops = ECEIDOps(SA[1.0 0.0; 0.0 0.0],sbm)
    RunECEID!(sbm,ops,NSTEPS,DT,store!, storage)
    return storage.sz
end

# Figure 1

sz = my_run(SBModel(OhmicSBParams(ϵ = 0.0, Δ = 1.0, ωc = 7.5, α = 0.1, No = 60, β = 5.0)))
refdata = readdlm("ECEID_ref_results/fig1_ECEID.dat")[:,2]
@test TOL > sum(abs.(refdata-sz))

# Figure 2

sz = my_run(SBModel(OhmicSBParams(ϵ = 0.0, Δ = 1.0, ωc = 7.5, α = 0.6, No = 60, β = 5.0)))
refdata = readdlm("ECEID_ref_results/fig2_ECEID.dat")[:,2]
@test TOL > sum(abs.(refdata-sz))

# Figure 3

sz = my_run(SBModel(OhmicSBParams(ϵ = 0.0, Δ = 1.0, ωc = 2.5, α = 1.2, No = 60, β = 0.2)))
refdata = readdlm("ECEID_ref_results/fig3_ECEID.dat")[:,2]
@test TOL > sum(abs.(refdata-sz))

# Figure 4

sz = my_run(SBModel(OhmicSBParams(ϵ = 1.0, Δ = 1.0, ωc = 7.5, α = 0.1, No = 60, β = 5.0)))
refdata = readdlm("ECEID_ref_results/fig4_ECEID.dat")[:,2]
@test TOL > sum(abs.(refdata-sz))

# Figure 5

sz = my_run(SBModel(OhmicSBParams(ϵ = 0.0, Δ = 1.0, ωc = 2.5, α = 0.09, No = 60, β = 0.1)))
refdata = readdlm("ECEID_ref_results/fig5_ECEID.dat")[:,2]
@test TOL > sum(abs.(refdata-sz))

# Figure 6

sz = my_run(SBModel(OhmicSBParams(ϵ = 0.0, Δ = 1.0, ωc = 2.5, α = 0.5625, No = 60, β = 10.0)))
refdata = readdlm("ECEID_ref_results/fig6_ECEID.dat")[:,2]
@test TOL > sum(abs.(refdata-sz))

# Figure 7

sz = my_run(SBModel(OhmicSBParams(ϵ = 0.0, Δ = 1.0, ωc = 2.5, α = 0.09, No = 60, β = 5.0)))
refdata = readdlm("ECEID_ref_results/fig7_ECEID.dat")[:,2]
@test TOL > sum(abs.(refdata-sz))


text:
using RCall
R"""
# install 'ggalt' pkg
# devtools::install_github("hrbrmstr/ggalt")
options(scipen = 999)
library(ggplot2)
library(ggalt)
midwest_select <- midwest[midwest$poptotal > 350000 &
                            midwest$poptotal <= 500000 &
                            midwest$area > 0.01 &
                            midwest$area < 0.1, ]

# Plot
ggplot(midwest, aes(x=area, y=poptotal)) +
  geom_point(aes(col=state, size=popdensity)) +   # draw points
  geom_smooth(method="loess", se=F) +
  xlim(c(0, 0.1)) +
  ylim(c(0, 500000)) +   # draw smoothing line
  geom_encircle(aes(x=area, y=poptotal),
                data=midwest_select,
                color="red",
                size=2,
                expand=0.08) +   # encircle
  labs(subtitle="Area Vs Population",
       y="Population",
       x="Area",
       title="Scatterplot + Encircle",
       caption="Source: midwest")
"""


text:
using BenchmarkTools
using Distributions
using GLM
using HypothesisTests
using LinearAlgebra
using Plots
using Random
using RDatasets
using SpecialFunctions
using Statistics

# Ridge regression

n = 10000
m = 1000

Random.seed!(666)
X = randn(n, m)

y = 10*X[:,1] + X[:,2] + randn(n)

# Exercise

ridge_reg(X, y, μ) = (X'*X + μ*I) \ (X'*y)

eigen_dec = eigen(X'*X)
Q = eigen_dec.vectors
Q_inv = Matrix(Q')
λ = eigen_dec.values

ridge_reg(X, y, μ, Q, Q_inv, λ) = Q * ((Diagonal(1 ./ (λ .+ μ)) * ( Q_inv * (X'*y))))

w1 = ridge_reg(X, y, 10)
w2 = ridge_reg(X, y, 10, Q, Q_inv, λ)

norm(w1 - w2)

# Comparison

@btime ridge_reg(X, y, 10);

@btime ridge_reg(X, y, 10, Q, Q_inv, λ);

# Dependence on mu

μs = range(0, 1000; length=50)

ws = hcat(ridge_reg.(Ref(X), Ref(y), μs, Ref(Q), Ref(Q_inv), Ref(λ))...)

plot(μs, abs.(ws');
    label="",
    yscale=:log10,
    xlabel="mu",
    ylabel="weights: log scale",
)

# LASSO

S(x, η) = max(x-η, 0) - max(-x-η, 0)

function lasso(X, y, μ, Q, Q_inv, λ;
        max_iter = 100,
        ρ = 1e3,
        w = zeros(size(X,2)),
        u = zeros(size(X,2)),
        z = zeros(size(X,2)),
    )

    for i in 1:max_iter
        w = Q * ( (Diagonal(1 ./ (λ .+ ρ)) * ( Q_inv * (X'*y + ρ*(z-u)))))
        z = S.(w + u, μ / ρ)
        u = u + w - z
    end
    return w, u, z
end

ws = zeros(size(X,2), length(μs))

for (i, μ) in enumerate(μs)
    global w, u, z
    w, u, z = i > 1 ? lasso(X, y, μ, Q, Q_inv, λ; w, u, z) : lasso(X, y, μ, Q, Q_inv, λ)
    ws[:,i] = w
end

plot(μs, abs.(ws');
    label="",
    yscale=:log10,
    xlabel="mu",
    ylabel="weights: log scale",
)

plot(0:0.1:10, gamma;
    xlabel="x",
    ylabel="gamma(x): log scale",
    label="",
    yscale=:log10,
)

# Exercise

volume_true(m, R) = π^(m/2) *R^2 / gamma(m/2 + 1)

plot(1:100, m -> volume_true.(m, 1);
    xlabel="dimension",
    ylabel="unit ball volume: log scale",
    label="",
    yscale=:log10,
)

# Exercise

function volume_monte_carlo(m::Int; n::Int=10000)
    X = 2*rand(m, n).-1
    X_norm_sq = sum(X.^2; dims=1)
    return 2^m*mean(X_norm_sq .<= 1)
end

# Volume of a ball

ms = 1:15
ns = Int64.([1e1; 1e3; 1e5])

Random.seed!(666)

plt = plot(ms, m -> volume_true(m, 1);
    xlabel="dimension",
    ylabel="unit ball volume",
    legend=:topleft,
    label="True",
    line=(4,:black),
)

for n in ns
    plot!(plt, ms, m -> volume_monte_carlo.(m; n=n); label="n = $n")
end

display(plt)

# Generating from the uniform distributions

rand(Uniform(-1, 1), 10, 5)

# Sampling from distributions

d1 = Normal()
d2 = Normal(1, 1)
d3 = Normal(0, 0.01)

f1(x) = pdf(d1, x)
f2(x) = pdf(d2, x)
f3(x) = pdf(d3, x)

function plot_histogram(xs, f; kwargs...)
    plt = histogram(xs;
        label="Sampled density",
        xlabel = "x",
        ylabel = "pdf(x)",
        nbins = 85,
        normalize = :pdf,
        opacity = 0.5,
        kwargs...
    )

    plot!(plt, range(minimum(xs), maximum(xs); length=100), f;
        label="True density",
        line=(4,:black),
    )

    return plt
end

plot_histogram(rand(d1, 1000000), f1)

# Exercise

function rejection_sampling(f, f_max, x_min, x_max; n=1000000)
    xs = x_min .+ (x_max - x_min)*rand(n)
    ps = f_max*rand(n)
    return xs[f.(xs) .>= ps]
end

# Rejection sampling

xlims = (-10, 10)

for (f, d) in zip((f1, f2, f3), (d1, d2, d3))
    Random.seed!(666)
    xs = rejection_sampling(f, f(d.μ), xlims...)

    pl = plot_histogram(xs, f)
    display(pl)
end

# Exercise

n = 1000
m = 9

Random.seed!(666)
xs = rand(m, n)

dist1 = [norm(x-y) for x in eachcol(xs), y in eachcol(xs)]
dist2 = [dist1[i,j] for i in 1:n for j in i+1:n]

extrema(dist2)

# Hypothesis testing

Random.seed!(666)

n = 1000
xs = randn(n)

# Exercise

t = mean(xs) / std(xs) * sqrt(n)

prob = cdf(TDist(n-1), t)
p = 2*min(prob, 1-prob)

# One-sample t-test

OneSampleTTest(xs)

# Linear models

wages = dataset("plm", "Snmesp")
wages.W = 2. .^ (wages.W)

X = Matrix(wages[:, [:N, :Y, :I, :K, :F]])
X = hcat(ones(size(X,1)), X)
y = wages[:, :W]

w0 = (X'*X) \ (X'*y)

model = lm(@formula(W ~ 1 + N + Y + I + K + F), wages)

# Exercise

norm(coef(model) - w0)

model_red = lm(@formula(W ~ 1 + N + Y + I + F), wages)

(r2(model), r2(model_red))

# Plot histograms

y_hat = predict(model)

plot_histogram(y_hat, x -> pdf(Normal(mean(y_hat), std(y_hat)), x))

plot_histogram(y_hat, x -> pdf(fit(Normal, y_hat), x))

test_normality = ExactOneSampleKSTest(y_hat, Normal(mean(y_hat), std(y_hat)))

# Generalized linear models

model = glm(@formula(W ~ 1 + N + Y + I + K + F), wages, InverseGaussian(), SqrtLink())

# Exercise

g_inv(z) = z^2

y_hat = g_inv.(X*coef(model))

scatter(y, y_hat;
    label="",
    xlabel="Label",
    ylabel="Prediction",
)

max_stars_repo_name: petrmvala/JuliaCourse

max_issues_repo_name: petrmvala/JuliaCourse

max_forks_repo_name: petrmvala/JuliaCourse


text:
# This file is a part of Julia. License is MIT: https://julialang.org/license

# RUN: julia --startup-file=no %s %t -O
# RUN: cat %t/* | FileCheck %s

include(joinpath("..", "testhelpers", "llvmpasses.jl"))

# JuliaLang/julia#38922
function haszerolayout(x::NTuple{32, VecElement{UInt8}})
    rx = Ref(x)
    GC.@preserve rx begin
        lower = iszero(unsafe_load(Ptr{UInt128}(pointer_from_objref(rx)), 1))
        upper = iszero(unsafe_load(Ptr{UInt128}(pointer_from_objref(rx)), 2))
        lower & upper
    end
end

# CHECK-LABEL: @julia_haszerolayout
# CHECK: top:
# CHECK-NOT: @jl_gc_pool_alloc
# CHECK: extractelement
# CHECK: ret i8
emit(haszerolayout, NTuple{32,VecElement{UInt8}})

max_stars_repo_name: jonas-schulze/julia

max_issues_repo_name: jonas-schulze/julia

max_forks_repo_name: jonas-schulze/julia


text:
# This file is auto-generated by AWSMetadata.jl
using AWS
using AWS.AWSServices: amp
using AWS.Compat
using AWS.UUIDs

"""
    create_workspace()
    create_workspace(params::Dict{String,<:Any})

Creates a new AMP workspace.

# Optional Parameters
Optional parameters can be passed as a `params::Dict{String,<:Any}`. Valid keys are:
- `"alias"`: An optional user-assigned alias for this workspace. This alias is for user
  reference and does not need to be unique.
- `"clientToken"`: Optional, unique, case-sensitive, user-provided identifier to ensure the
  idempotency of the request.
"""
create_workspace(; aws_config::AbstractAWSConfig=global_aws_config()) = amp("POST", "/workspaces", Dict{String, Any}("clientToken"=>string(uuid4())); aws_config=aws_config)
create_workspace(params::AbstractDict{String, <:Any}; aws_config::AbstractAWSConfig=global_aws_config()) = amp("POST", "/workspaces", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("clientToken"=>string(uuid4())), params)); aws_config=aws_config)

"""
    delete_workspace(workspace_id)
    delete_workspace(workspace_id, params::Dict{String,<:Any})

Deletes an AMP workspace.

# Arguments
- `workspace_id`: The ID of the workspace to delete.

# Optional Parameters
Optional parameters can be passed as a `params::Dict{String,<:Any}`. Valid keys are:
- `"clientToken"`: Optional, unique, case-sensitive, user-provided identifier to ensure the
  idempotency of the request.
"""
delete_workspace(workspaceId; aws_config::AbstractAWSConfig=global_aws_config()) = amp("DELETE", "/workspaces/$(workspaceId)", Dict{String, Any}("clientToken"=>string(uuid4())); aws_config=aws_config)
delete_workspace(workspaceId, params::AbstractDict{String, <:Any}; aws_config::AbstractAWSConfig=global_aws_config()) = amp("DELETE", "/workspaces/$(workspaceId)", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("clientToken"=>string(uuid4())), params)); aws_config=aws_config)

"""
    describe_workspace(workspace_id)
    describe_workspace(workspace_id, params::Dict{String,<:Any})

Describes an existing AMP workspace.

# Arguments
- `workspace_id`: The ID of the workspace to describe.

"""
describe_workspace(workspaceId; aws_config::AbstractAWSConfig=global_aws_config()) = amp("GET", "/workspaces/$(workspaceId)"; aws_config=aws_config)
describe_workspace(workspaceId, params::AbstractDict{String, <:Any}; aws_config::AbstractAWSConfig=global_aws_config()) = amp("GET", "/workspaces/$(workspaceId)", params; aws_config=aws_config)

"""
    list_workspaces()
    list_workspaces(params::Dict{String,<:Any})

Lists all AMP workspaces, including workspaces being created or deleted.

# Optional Parameters
Optional parameters can be passed as a `params::Dict{String,<:Any}`. Valid keys are:
- `"alias"`: Optional filter for workspace alias. Only the workspaces with aliases that
  begin with this value will be returned.
- `"maxResults"`: Maximum results to return in response (default=100, maximum=1000).
- `"nextToken"`: Pagination token to request the next page in a paginated list. This token
  is obtained from the output of the previous ListWorkspaces request.
"""
list_workspaces(; aws_config::AbstractAWSConfig=global_aws_config()) = amp("GET", "/workspaces"; aws_config=aws_config)
list_workspaces(params::AbstractDict{String, Any}; aws_config::AbstractAWSConfig=global_aws_config()) = amp("GET", "/workspaces", params; aws_config=aws_config)

"""
    update_workspace_alias(workspace_id)
    update_workspace_alias(workspace_id, params::Dict{String,<:Any})

Updates an AMP workspace alias.

# Arguments
- `workspace_id`: The ID of the workspace being updated.

# Optional Parameters
Optional parameters can be passed as a `params::Dict{String,<:Any}`. Valid keys are:
- `"alias"`: The new alias of the workspace.
- `"clientToken"`: Optional, unique, case-sensitive, user-provided identifier to ensure the
  idempotency of the request.
"""
update_workspace_alias(workspaceId; aws_config::AbstractAWSConfig=global_aws_config()) = amp("POST", "/workspaces/$(workspaceId)/alias", Dict{String, Any}("clientToken"=>string(uuid4())); aws_config=aws_config)
update_workspace_alias(workspaceId, params::AbstractDict{String, <:Any}; aws_config::AbstractAWSConfig=global_aws_config()) = amp("POST", "/workspaces/$(workspaceId)/alias", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("clientToken"=>string(uuid4())), params)); aws_config=aws_config)