text:
# ------------------------------------------------------------------
# Licensed under the MIT License. See LICENSE in the project root.
# ------------------------------------------------------------------

"""
    intersecttype(f, b1, b2)

Compute the intersection type of two boxes `b1` and `b2`
and apply function `f` to it.

The intersection type can be one of four types:

1. overlap with non-zero measure
2. intersect at one of the boundaries
3. intersect at corner point
4. do not overlap nor intersect
"""
function intersecttype(f::Function, b1::Box{Dim,T}, b2::Box{Dim,T}) where {Dim,T}
  m1, M1 = coordinates.(extrema(b1))
  m2, M2 = coordinates.(extrema(b2))

  # relevant vertices
  u = Point(max.(m1, m2))
  v = Point(min.(M1, M2))

  # branch on possible configurations
  if u ≺ v
    return OverlappingBoxes(Box(u, v)) |> f
  elseif u ≻ v
    return NoIntersection() |> f
  elseif isapprox(u, v, atol=atol(T))
    return CornerTouchingBoxes(u) |> f
  else
    return FaceTouchingBoxes(Box(u, v)) |> f
  end
end

"""
    intersecttype(function, ray, box)

Calculate the intersection type of a `ray`` and a `box`` and apply function `f` to it.

## References

* Williams A, Barrus S, Morley R K, et al., 2005. [An efficient and robust ray-box
  intersection algorithm]
  (https://dl.acm.org/doi/abs/10.1145/1198555.1198748)
"""
function intersecttype(f::Function, r::Ray{Dim,T}, b::Box{Dim,T}) where {Dim,T}
  invdir = one(T) ./ direction(r)
  lo, up = coordinates.(extrema(b))
  orig = coordinates(origin(r))

  tmin = zero(T)
  tmax = typemax(T)

  # check for intersection with slabs along with each axis
  for i in 1:Dim
    imin = (lo[i] - orig[i]) * invdir[i]
    imax = (up[i] - orig[i]) * invdir[i]

    # swap variables if necessary
    invdir[i] < zero(T) && ((imin, imax) = (imax, imin))

    # the ray is on a face of the box, avoid NaN
    (isnan(imin) || isnan(imax)) && continue

    (tmin > imax || imin > tmax) && return NoIntersection() |> f

    tmin = max(tmin, imin)
    tmax = min(tmax, imax)
  end

  tmin ≈ tmax && return TouchingRayBox(r(tmin)) |> f

  return CrossingRayBox(Segment(r(tmin), r(tmax))) |> f
end


text:
println("Running tests for Jags-j0.5-v1.0.2:")

using Jags
using Base.Test

old = pwd()

using Base.Test

code_tests = ["test_cmd.jl";]

println("Run execution_tests only if Jags is available")
execution_tests = [
  "test_line1.jl",
  "test_line2.jl",
  "test_rats.jl",
  "test_bones1.jl",
  "test_bones2.jl",
  "test_dyes.jl";
]

for my_test in code_tests
    println("\n  * $(my_test) *")
    include(my_test)
end

if isdefined(Main, :JAGS_HOME) && length(JAGS_HOME) > 0
  try
    for my_test in execution_tests
        println("\n  * $(my_test) *")
        include(my_test)
    end
    println()
  catch e
     println("Is Jags properly installed?")
     println(e)
     println("No simulation runs have been performed.")
  end
else
  println("\n\nJAGS_HOME not found. Skipping all tests that depend on the Jags executable!\n")
end

cd(old)

ext: jl


text:
"""
    getasymptotes(approx::AbstractMMAApprox)

Returns the asymptotes of the MMA approximation `approx`.
"""
function getasymptotes(approx::AbstractMMAApprox)
    throw("`getasymptotes` is not defined.")
end

"""
    getapproxfg(approx::AbstractMMAApprox)

Returns the last cached value of the output of `approx`.
"""
function getapproxfg(approx::AbstractMMAApprox)
    throw("`getapproxfg` is not defined.")
end

docσ = """
    getσ(approx::AbstractMMAApprox)
    setσ!(approx::AbstractMMAApprox, σ)

Get or set the move limit `σ` of the MMA approximation `approx`. See [`MMAApprox`](@ref) for an explanation.
"""

@doc docσ
function getσ(approx::AbstractMMAApprox)
    throw("`getσ` is not defined.")
end

@doc docσ
function setσ!(approx::AbstractMMAApprox, σ)
    throw("`setσ!` is not defined.")
end

docρ = """
    getρ(approx::AbstractMMAApprox)
    setρ!(approx::AbstractMMAApprox, ρ)

Get or set the value of `ρ` in `approx`. See [`MMAApprox`](@ref) for an explanation.
"""

@doc docρ
function getρ(approx::AbstractMMAApprox)
    throw("`getρ` is not defined.")
end

@doc docρ
function setρ!(approx::AbstractMMAApprox, ρ)
    throw("`setρ!` is not defined.")
end

docfk = """
    getfk(approx::AbstractMMAApprox)
    setfk!(approx::AbstractMMAApprox, fk)

Get or set the function value `fk` of the approximated function at the approximation point `xk` of `approx`.
"""

@doc docfk
function getfk(approx::AbstractMMAApprox)
    throw("`getfk` is not defined.")
end

@doc docfk
function setfk!(approx::AbstractMMAApprox, f)
    throw("`setfk!` is not defined.")
end

docxk = """
    getxk(approx::AbstractMMAApprox)
    setxk!(approx::AbstractMMAApprox, x::AbstractVector)

Get or set the approximation point `xk` of `approx`.
"""

@doc docxk
function getxk(approx::AbstractMMAApprox)
    throw("`getfk` is not defined.")
end

@doc docxk
function setxk!(approx::AbstractMMAApprox, x::AbstractVector)
    throw("`setxk!` is not defined.")
end

doc∇fk = """
    get∇fk(approx)
    set∇fk!(approx, ∇f::AbstractVecOrMat)

Get or set the value of the gradient/Jacobian, `∇fk`, of the approximated function at the approximation point `xk`.
"""

@doc doc∇fk
function get∇fk(approx::AbstractMMAApprox)
    throw("`get∇fk` is not defined.")
end

@doc doc∇fk
function set∇fk!(approx::AbstractMMAApprox, ∇f::AbstractVecOrMat)
    throw("`set∇fk!` is not defined.")
end

docupdateapprox! = """
    updateapprox!(f, x)

Updates the MMA approximation in `f` making it around the point `x`.
"""

@doc docupdateapprox!
function updateapprox!(approx::AbstractMMAApprox, x::AbstractVector)
    throw("`updateapprox` is not defined.")
end

"""
    updateapprox!(approx::AbstractMMAApprox, x, f, ∇f)

Updates the MMA approximation `approx` making it around the point `x` where `f` is the value of the approximated function at `x` and `∇f` is its gradient if `f` is a number and the Jacobian if `f` is a vector.
"""
function updateapprox!(approx::AbstractMMAApprox, x, f, ∇f)
    throw("`updateapprox` is not defined.")
end

"""
    updateapprox!(approx::AbstractMMAApprox)

Updates the MMA approximation `approx`. This is called when `σ`, `ρ` or `xk` change.
"""
function updateapprox!(approx::AbstractMMAApprox)
    throw("`updateapprox` is not defined.")
end


text:
using Pkg

Pkg.develop(path = "..")
# this is needed since Publish v0.9 breaks our theming hack
Pkg.pin(name = "Publish", version = "0.8")

using Publish
using Pkg.Artifacts
using Metalhead

# override default theme
Publish.Themes.default() = artifact"flux-theme"

p = Publish.Project(Metalhead)

# needed to prevent error when overwriting
rm("dev", recursive = true, force = true)
rm(p.env["version"], recursive = true, force = true)

# build documentation
deploy(Metalhead; root = "/Metalhead.jl", force = true, label = "dev")


text:
# used to generate error codes in src/error_codes.jl
# if run with no args, prints the generated code to stdout
# with a file argument, write the code to that file
# e.g., julia deps/error_codes.jl src/error_codes.jl
#
# only needs to be re-run if the Travis diff job fails

using EzXML
using HTTP

const ERROR_CODE_APPENDIX = "https://www.postgresql.org/docs/current/errcodes-appendix.html"
const HEADER = "# autogenerated by deps/error_codes.jl using $ERROR_CODE_APPENDIX"

const SUCCESS_CLASS = "C00"
const WARNING_CLASSES = ("C01", "C02")
const EXTERNAL_CLASSES = ("C38", "C39")

pascalcase(str) = replace(titlecase(str), '_' => "")

error_code_html(url=ERROR_CODE_APPENDIX) = String(HTTP.get(url))

function error_code_table(html)
    parsed = parsehtml(html)

    return findfirst("//table[@summary='PostgreSQL Error Codes']", parsed)
end

parse_row(node) = map(nodecontent, findall("td/code", node))

function generate_error_codes(io, html=error_code_html())
    table = error_code_table(html)
    rows = findall("tbody/tr/td/code/../..", table)
    id_names = Set{String}()

    class_enum_io = IOBuffer()
    println(class_enum_io, "@cenum(\n    # Enum type\n    Class,\n    # Enum values")
    error_enum_io = IOBuffer()
    println(error_enum_io, "@cenum(\n    # Enum type\n    ErrorCode,\n    # Enum values")
    alias_io = IOBuffer()
    error_names_io = IOBuffer()
    println(error_names_io, "const ERROR_NAMES = Dict(")

    for row in rows
        code, name = parse_row(row)

        id_name = pascalcase(name)
        class = "C$(code[1:2])"
        error_code = "E$code"

        if endswith(code, "000")
            suffix = class in (SUCCESS_CLASS, WARNING_CLASSES...) ? "Class" : "ErrorClass"

            println(class_enum_io, "    $class,")
            println(
                alias_io,
                "\n# $id_name\nconst $(id_name)$(suffix) = PQResultError{$class}\n",
            )
        end

        if class in WARNING_CLASSES
            id_name *= "Warning"
        end

        if id_name in id_names
            if class in EXTERNAL_CLASSES
                id_name *= "Ext"  # like psycopg2
            else
                error("Duplicate identifier: $id_name")
            end
        end

        push!(id_names, id_name)

        println(alias_io, "const $(id_name) = PQResultError{$class,$error_code}")
        println(error_enum_io, "    $error_code,")
        println(error_names_io, "    $id_name => \"$id_name\",")
    end

    # unknown error
    class = "CUN"
    error_code = "EUNOWN"
    println(class_enum_io, "    $class,")
    println(
        alias_io, "\n# Unknown Error\nconst UnknownErrorClass = PQResultError{$class}\n"
    )
    println(alias_io, "const UnknownError = PQResultError{$class,$error_code}")
    println(error_enum_io, "    $error_code,")
    println(error_names_io, "    UnknownError => \"UnknownError\",")

    println(class_enum_io, ")")
    println(error_enum_io, ")")
    println(error_names_io, ")")

    foreach(seekstart, (class_enum_io, error_enum_io, alias_io, error_names_io))
    write(
        io,
        HEADER,
        "\n\n",
        class_enum_io,
        "\n",
        error_enum_io,
        alias_io,
        "\n",
        error_names_io,
    )

    return
end

function main()
    if length(ARGS) > 0
        open(ARGS[1], "w") do fp
            generate_error_codes(fp)
        end
    else
        generate_error_codes(stdout)
    end
end

abspath(PROGRAM_FILE) == @__FILE__() && main()


text:
# Author: Hiroharu Sugawara <hsugawa@gmail.com>
# Copyright: (C) 2020 Hiroharu Sugawara

export investigate_image, investigate_codeblock

function newmodule(name::String, outdir::String)::Module
   mod = Core.eval(Main, Meta.parse("""
module $name
   macro OUTPUT(); return "$outdir"; end
end
"""))
end

function run_code(mod, code, out_path, err_path)
   function run_code_inside(mod, code, out_path, err_path)
      res=nothing
      open(out_path, "w") do out_f
         redirect_stdout(out_f) do
            open(err_path, "w") do err_f
               n = sizeof(code)
               pos = 1
               while pos <= n
                  try
                     ex, pos = Meta.parse(code, pos)
                     isnothing(ex) && continue
                     res = Core.eval(mod, ex)
                  catch exc
                     if isa(exc,Base.Meta.ParseError)
                        print(err_f, "ERROR: syntax: "*exc.msg)
                     else
                        showerror(err_f, exc)
                     end
                     break
                  end
               end
            end # open(err_path)
         end # redirect_stdout(outf)
      end # open(out_path
      return res
   end

   isempty(code) && return nothing
   ispath(out_path) || mkpath(dirname(out_path))
   ispath(err_path) || mkpath(dirname(err_path))
   run_code_inside(mod,code, out_path, err_path)
end

function remove_hide_from_code(code)
   removed=String[]
   for text in split(code,"\n")
      endswith(text,"#hide") && continue
      push!(removed,text)
   end
   removed=join(removed,"\n")
   isempty(removed) || return removed
   endswith(removed,"\n") && (removed *= "\n")
   return removed
end

OUTDIR="assets"
global MOD=newmodule("mod", OUTDIR)
global COUNT=0

function run_the_code(code)
   global COUNT, OUTDIR, MOD
   COUNT += 1
   out_path=joinpath(OUTDIR, "code"*string(COUNT)*".out")
   err_path=joinpath(OUTDIR, "code"*string(COUNT)*".err")
   res=run_code(MOD, code, out_path, err_path)
   out_message=join(readlines(out_path),"\n")
   err_message=join(readlines(err_path),"\n")
   return res, out_message, err_message
end

function destruct_attrs(a)
   id=""
   length(a) >= 1 && (id=a[1])
   classes=[]
   length(a) >= 2 && (classes=a[2])
   keyvals=[]
   length(a) >= 3 && (keyvals=a[3])
   return id, classes, keyvals
end

function investigate_codeblock(t, c, format, meta)
   function stringfy(x)
      # isa(x, AbstractString) && return(x)
      io=IOBuffer()
      show(io, "text/plain", x)
      String(take!(io))
   end

   (t == "CodeBlock") || return nothing
   attrs, code=c
   id, classes, keyvals=destruct_attrs(attrs)
   length(classes)>=1 || return nothing
   class1 = classes[1]
   startswith(class1,"julia-exec,") || return nothing
   class1 = replace(class1, "julia-exec," => "")
   res, out_s, err_s = run_the_code(code)
   attrs1=deepcopy(attrs)
   attrs1[2][1]=class1
   code=remove_hide_from_code(code)
   isempty(code) && return []
   element1=CodeBlock(attrs1,code)
   result=Any[]
   push!(result, element1)

   if ! isempty(out_s)
      element2=CodeBlock(attrs1,out_s)
      push!(result, element2)
   end
   if ! isempty(err_s)
      element3=CodeBlock(attrs1,err_s)
      push!(result, element3)
   end
   #
   if isempty(out_s) && isempty(err_s)
      if !isnothing(res)
         element4=CodeBlock(attrs1,stringfy(res))
         push!(result, element4)
      end
   end
   return result
end

function investigate_image(t, c, format, meta)
   (t == "Image") || return nothing
   attr, inlines, url_title = c
   url, title = url_title
   url = replace(url, "@OUTPUT/" => "assets/")
   Image(attr, inlines, [url, title])
end


text:
using Documenter, Oracle

makedocs(
    sitename="Oracle.jl",
    pages = [ "Home" => "index.md",
              "Tutorial" => "tutorial.md",
              "API Reference" => "api.md"
            ],
)

deploydocs(
    repo = "github.com/felipenoris/Oracle.jl.git",
    target = "build",
)


text:
using Plots
include("RRT.jl")
closeall()
plotly()

cost, isPathFound, nlist = rrtPathPlanner(30)
@show isPathFound
@show cost
@show length(nlist)
@show nlist

iterList = [40 80 100 120 150 180 200 300]
nTrials = 30
cost = 0
listCosts = Vector{Float32}()
listpSucc= Vector{Float32}()

print("\n\n")
for testIter in iterList
    totalcost =0
    idx = 0
    @show testIter
    nSuccess = 0
    cost = 0
    while idx < nTrials
        cost, isPathFound, nlist = rrtPathPlanner(testIter)
        if isPathFound
            nSuccess += 1
        end
        totalcost += cost
        idx += 1
    end
    avgCost = totalcost /  nSuccess
    @show nSuccess
    pSucc = nSuccess / nTrials
    @show avgCost
    @show listCosts
    @show listpSucc
    print("\n\n")
    #@printf("For the iter of %d the avg cost was %d across %d trials", maxIter, avgCost, nTrials)
    push!(listCosts, avgCost)
    push!(listpSucc, pSucc)
end

@show iterList
@show listCosts
@show listpSucc

sizeplot = (400, 600)

# plot nIter vs pathcost
 pRRTcost = plot(iterList, listCosts', show=true, size=sizeplot, seriestype=:scatter,
      legend=false, yaxis=((0,100), 0:20:100), xaxis=((0,320), 0:50:300), color=:black)

title!("RRT nIter vs pathcost")
ylabel!("euclidean path cost")

# plot nIter vs success rate
pRRTsuccess = scatter(iterList, listpSucc', legend=false, size=sizeplot, yaxis=((0,1), 0:0.1:1), xaxis=((0,320), 0:50:300), color=:orange)
title!("RRT nIter vs pSuccess")
ylabel!("# P(success) = numSucc / 30 trials")
xlabel!("iterations requested of RRT alg")

plot(pRRTcost, pRRTsuccess, layout=(2,1),legend=false)

# [slack] <mkborregaard> twinx
# [slack] <mkborregaard> plot(something); plot!(twinx(), something)

########################################
### DO NOT USE (why is this here?)
#pRRTsuccess[:axis]

# y = rand(100)
# plot(0:10:100,rand(11,4),lab="lines",w=3,palette=:grays,fill=0,α=0.6)
# scatter!(y,zcolor=abs(y - 0.5),m=(:heat,0.8,stroke(1,:green)),ms=10 * abs(y - 0.5) + 4,lab="grad")
 #
 #
# cost, isPathFound, nlist = rrtPathPlanner(40) #maxIter
# plotPath(isPathFound,nlist)
 #
#


text:
using ForwardDiff
using Parameters: @with_kw, @unpack
using Setfield: set
using StaticArrays: SVector

import ..Continuations: get_prob_cache, get_u0, residual!, residual_jacobian!,
    isindomain

using ..Codim1LimitCycle
using ..Codim1LimitCycle: LimitCycleProblem, LimitCycleCache, residual_lc!

abstract type
    Codim2LCProblem{skind, tkind} <: BifurcationProblem{skind, tkind}
end

struct FoldLimitCycleProblem{
        skind <: MutableState,  # <: StateKind,
        tkind <: TimeKind,
        P <: LimitCycleProblem,
        XS <: AbstractMatrix,
        R <: Real,
        T, L1, L2,
        } <: Codim2LCProblem{skind, tkind}
    statekind::skind
    timekind::tkind
    super::P
    xs0::XS
    vs0::XS
    l0::R
    dl0::R
    t0::T
    t_domain::Tuple{T, T}
    param_axis1::L1
    param_axis2::L2
end

function FoldLimitCycleProblem(
        prob::LimitCycleProblem;
        vs0,
        dl0,
        xs0 = prob.xs0,
        l0 = prob.l0,
        param_axis1 = prob.param_axis,
        param_axis2,
        t0 = SVector(get(prob.de_prob.p, param_axis1),
                     get(prob.de_prob.p, param_axis2)),
        t2_domain = (typemin(eltype(t0)), typemax(eltype(t0))),
        t_domain = (
            SVector(prob.t_domain[1], t2_domain[1]),
            SVector(prob.t_domain[2], t2_domain[2]),
        ),
        )
    return FoldLimitCycleProblem(
        MutableState(),
        timekind(prob.de_prob),
        prob,
        xs0,
        vs0,
        l0,
        dl0,
        t0,
        t_domain,
        param_axis1,
        param_axis2,
    )
end

function FoldLimitCycleProblem(;
        vs0,
        t0,
        param_axis1,
        param_axis2,
        kwargs...)
    prob = LimitCycleProblem(;
        t0 = t0[1],
        param_axis = param_axis1,
        kwargs...)
    return FoldLimitCycleProblem(
        prob;
        vs0 = vs0,
        t0 = t0,
        param_axis1 = param_axis1,
        param_axis2 = param_axis2,
    )
end

# ----------------------------------------------------------------------- cache

struct FoldLimitCycleCache{
        P <: FoldLimitCycleProblem,
        C <: LimitCycleCache,
        } <: AbstractProblemCache{P}
    super::C
    prob::P
end

function FoldLimitCycleCache(prob)
    return FoldLimitCycleCache(
        LimitCycleCache(prob.super),
        prob,
    )
end

# ------------------------------------------------------ continuation interface

get_prob_cache(prob::FoldLimitCycleProblem) = FoldLimitCycleCache(prob)

get_u0(prob::FoldLimitCycleProblem) = vcat(
    view(prob.xs0, :),          # \__ ξ
    prob.l0,                    # /
    view(prob.vs0, :),          # \__ η
    prob.dl0,                   # /
    prob.t0,                    # --- α
)
# modified_param! requires prob.t0 to be at the end

Codim1LimitCycle.u_idx_param(cache::FoldLimitCycleCache) =
    2(length(cache.prob.xs0) + 1) .+ (1:2)

isindomain(u, cache::FoldLimitCycleCache) =
    Codim1LimitCycle.isindomain_lc(u, cache)

# ------------------------------------------------------------------- residual!

# H = [
#     H_lc = [                                           # F(ξ, α) = 0
#         state_condition
#         phase_condition
#     ]
#     eigvec_condition                                   # ∂₁ F(ξ, α) η = 0
#     eigvec_constraint                                  # <η,η> - 1 = 0
# ]

# Almost identical to: [[../codim2/diffeq.jl::modified_param!]]
function modified_param!(p, u)
    t1, t2 = u[end-1:end]
    q0 = p.super.de_prob.p
    q1 = set(q0, p.param_axis1, t1)
    q2 = set(q1, p.param_axis2, t2)
    return q2
end

function residual!(H, u, cache::FoldLimitCycleCache)
    prob = cache.prob
    q = modified_param!(prob, u)

    n_lc = length(prob.xs0) + 1
    ξ = view(u, 1:n_lc)                  # [xs; period]
    η = view(u, (1:n_lc) .+ n_lc)        # ∂ξ
    H_lc = view(H, 1:n_lc)               # F(ξ, α)
    H_ev = view(H, (1:n_lc) .+ n_lc)     # ∂₁ F(ξ, α) η
    H_ec = @view H[end]                  # <η,η> - 1

    ForwardDiff.jacobian!(
        reshape(H_ev, :, 1),  # =  ∂₁ F(ξ, α) η  =  ∂ᵧ F(ξ + γ η, α)
        (H_lc, d) -> residual_lc!(H_lc, (@. ξ + d[1] * η), q, cache.super),
        H_lc,
        SVector(zero(eltype(η))),
        # TODO: setup cache
    )

    H_ec[] = η ⋅ η - 1

    return H
end

# ---------------------------------------------------------- residual_jacobian!

# J = [
#     ∂₁F                              0       ∂₂F
#     (∂ᵪ(∂ᵧF(χ + γ η, α)))(χ = ξ)    ∂₁F    ∂₂((∂₁F)η)
#     0                                 2η      0
# ]

function residual_jacobian!(H, J, u, cache::FoldLimitCycleCache)
    prob = cache.prob
    q = modified_param!(prob, u)

    n_lc = length(prob.xs0) + 1
    ξ = view(u, 1:n_lc)                  # [xs; period]
    η = view(u, (1:n_lc) .+ n_lc)        # ∂ξ
    H_lc = view(H, 1:n_lc)               # F(ξ, α)
    H_ev = view(H, (1:n_lc) .+ n_lc)     # ∂₁ F(ξ, α) η

    fill!(J, zero(eltype(J)))

    @views H[end] = η ⋅ η - 1  # H_ec = <η,η> - 1
    @views J[end, n_lc .+ (1:n_lc)] .= 2 .* η

    ∂₁F = view(J, 1:n_lc, 1:n_lc)  # = ∂₁F(ξ, α)
    ForwardDiff.jacobian!(
        ∂₁F,
        (H_lc, ξ) -> residual_lc!(H_lc, ξ, q, cache.super),
        H_lc,
        ξ,
        # TODO: setup cache
    )
    @views J[n_lc .+ (1:n_lc), n_lc .+ (1:n_lc)] .= ∂₁F

    ForwardDiff.jacobian!(
        view(J, n_lc .+ (1:n_lc), 1:n_lc),  # (∂ᵪ(∂ᵧF(χ + γ η, α)))(χ = ξ)
        (y, x) -> ForwardDiff.jacobian!(
            reshape(view(y, :), :, 1),
            (H_lc, d) -> residual_lc!(H_lc, (@. x + d[1] * η), q, cache.super),
            fill!(similar(y), zero(eltype(y))),  # TODO: pre-allocate
            SVector(zero(eltype(η))),
        ),
        H_ev,  # y = ∂₁ F(ξ, α) η
        ξ,     # x
    )

    # Compute ∂₂F and ∂₂((∂₁F)η) in one go:
    ForwardDiff.jacobian!(
        (@view J[1:2n_lc, end-1:end]),  # [∂₂F; ∂₂((∂₁F)η)]
        (H_, α) -> begin
            p = modified_param!(prob, α)
            ForwardDiff.jacobian!(
                reshape(view(H_, (1:n_lc) .+ n_lc), :, 1),  # H_ev
                (H_lc, d) -> residual_lc!(H_lc, (@. ξ + d[1] * η),
                                          p, cache.super),
                view(H_, 1:n_lc),  # H_lc
                SVector(zero(eltype(η))),
            )
        end,
        view(H, 1:2n_lc),  # = H_ = [H_lc; H_ev]
        u[end-1:end],      # = α = parameters
        # TODO: setup cache
    )
    return (H, J)
end
# TODO: Some computations (e.g., H_lc) are done multiple times.  Don't.


text:
using ISICArchive
using Test
using JSON

# Basic tests that see if calls suceed.
# i.e. these test will detect drastic api changes

dataset_list = get(DatasetListRequest(limit = 4))
@test length(dataset_list) == 4

info = get(DatasetMetadataRequest(dataset_list[1]))
info2 = ISICArchive.DatasetMetadata(JSON.parse(json(info)))

@test info == info
@test info >= info
@test info <= info
@test (info != info) == false
@test (info > info) == false
@test (info < info) == false

@test info == info2
@test info >= info2
@test info <= info2
@test (info != info2) == false
@test (info > info2) == false
@test (info < info2) == false

image_list = get(ImageListRequest(datasetId = dataset_list[1], limit = 5))
@test length(image_list) == 5
image_list = get(ImageListRequest(info, limit = 5))
@test length(image_list) == 5
image_list = get(ImageListRequest(dataset_list[1], limit = 5))
@test length(image_list) == 5

img_info = get(ImageMetadataRequest(id = "5592ac579fc3c13155a57a80"))
img_info = get(ImageMetadataRequest(image_list[1]))
img_info2 = ISICArchive.ImageMetadata(JSON.parse(json(img_info)))

img = get(ImageThumbnailRequest(id = "5592ac579fc3c13155a57a80"))
img = get(ImageThumbnailRequest(img_info))
img = get(ImageThumbnailRequest(image_list[1]))

# download of high-quality will not be tested


text:
"""
    add_player!(scores, name, score=0)

Add a player to the score board.
"""
add_player!(scores, name, score=0) = push!(scores, name => score)

"""
    remove_player!(scores, name)

Remove a player from the score board.
"""
remove_player!(scores, name) = delete!(scores, name)

"""
    update_score!(scores, name, score)

Update the score of a player on the score board.
"""
update_score!(scores, name, score) = scores[name] = score

"""
    reset_score!(scores, name)

Reset the score of a player on the score board.
"""
reset_score!(scores, name) = update_score!(scores, name, 0)

"""
    players(scores)

Return a list of all players on the score board.
"""
players(scores) = collect(keys(scores))


text:
immutable SMatrix{S1, S2, T, L} <: StaticMatrix{T}
    data::NTuple{L, T}

    function SMatrix(d::NTuple{L,T})
        check_smatrix_params(Val{S1}, Val{S2}, T, Val{L})
        new(d)
    end

    function SMatrix(d::NTuple{L})
        check_smatrix_params(Val{S1}, Val{S2}, T, Val{L})
        new(convert_ntuple(T, d))
    end
end

@generated function check_smatrix_params{S1,S2,L}(::Type{Val{S1}}, ::Type{Val{S2}}, T, ::Type{Val{L}})
    if !(T <: DataType) # I think the way types are handled in generated fnctions might have changed in 0.5?
        return :(error("SMatrix: Parameter T must be a DataType. Got $T"))
    end

    if !isa(S1, Int) || !isa(S2, Int) || !isa(L, Int) || S1 < 0 || S2 < 0 || L < 0
        return :(error("SMatrix: Sizes must be positive integers. Got $S1 × $S2 ($L elements)"))
    end

    if S1*S2 == L
        return nothing
    else
        str = "Size mismatch in SMatrix. S1 = $S1, S2 = $S2, but recieved $L elements"
        return :(error($str))
    end
end

@generated function (::Type{SMatrix{S1}}){S1,L}(x::NTuple{L})
    S2 = div(L, S1)
    if S1*S2 != L
        error("Incorrect matrix sizes. $S1 does not divide $L elements")
    end
    T = promote_tuple_eltype(x)

    return quote
        $(Expr(:meta, :inline))
        SMatrix{S1, $S2, $T, L}(x)
    end
end

@generated function (::Type{SMatrix{S1,S2}}){S1,S2,L}(x::NTuple{L})
    T = promote_tuple_eltype(x)

    return quote
        $(Expr(:meta, :inline))
        SMatrix{S1, S2, $T, L}(x)
    end
end

@generated function (::Type{SMatrix{S1,S2,T}}){S1,S2,T,L}(x::NTuple{L})
    return quote
        $(Expr(:meta, :inline))
        SMatrix{S1, S2, T, L}(x)
    end
end

@inline convert{S1,S2,T}(::Type{SMatrix{S1,S2}}, a::StaticArray{T}) = SMatrix{S1,S2,T}(Tuple(a))
@inline convert{S1,S2,T}(::Type{SMatrix{S1,S2}}, a::AbstractArray{T}) = SMatrix{S1,S2,T}((a...))
@inline SMatrix(a::StaticMatrix) = SMatrix{size(typeof(a),1),size(typeof(a),2)}(Tuple(a))

#=
@inline (::Type{SMatrix{S1}}){S1}(x1) = SMatrix{S1}((x1,))
@inline (::Type{SMatrix{S1}}){S1}(x1,x2) = SMatrix{S1}((x1,x2))
@inline (::Type{SMatrix{S1}}){S1}(x1,x2,x3) = SMatrix{S1}((x1,x2,x3))
@inline (::Type{SMatrix{S1}}){S1}(x1,x2,x3,x4) = SMatrix{S1}((x1,x2,x3,x4))
@inline (::Type{SMatrix{S1}}){S1}(x...) = SMatrix{S1}(x)

@inline (::Type{SMatrix{S1,S2}}){S1,S2}(x1) = SMatrix{S1,S2}((x1,))
@inline (::Type{SMatrix{S1,S2}}){S1,S2}(x1,x2) = SMatrix{S1,S2}((x1,x2))
@inline (::Type{SMatrix{S1,S2}}){S1,S2}(x1,x2,x3) = SMatrix{S1,S2}((x1,x2,x3))
@inline (::Type{SMatrix{S1,S2}}){S1,S2}(x1,x2,x3,x4) = SMatrix{S1,S2}((x1,x2,x3,x4))
@inline (::Type{SMatrix{S1,S2}}){S1,S2}(x...) = SMatrix{S1,S2}(x)

@inline (::Type{SMatrix{S1,S2,T}}){S1,S2,T}(x1) = SMatrix{S1,S2,T}((x1,))
@inline (::Type{SMatrix{S1,S2,T}}){S1,S2,T}(x1,x2) = SMatrix{S1,S2,T}((x1,x2))
@inline (::Type{SMatrix{S1,S2,T}}){S1,S2,T}(x1,x2,x3) = SMatrix{S1,S2,T}((x1,x2,x3))
@inline (::Type{SMatrix{S1,S2,T}}){S1,S2,T}(x1,x2,x3,x4) = SMatrix{S1,S2,T}((x1,x2,x3,x4))
@inline (::Type{SMatrix{S1,S2,T}}){S1,S2,T}(x...) = SMatrix{S1,S2,T}(x)
=#

# Some more advanced constructor-like functions
@inline eye{N}(::Type{SMatrix{N}}) = eye(SMatrix{N,N})
@inline eye{N,M}(::Type{SMatrix{N,M}}) = eye(SMatrix{N,M,Float64})
@inline zeros{N,M}(::Type{SMatrix{N,M}}) = zeros(SMatrix{N,M,Float64})
@inline ones{N,M}(::Type{SMatrix{N,M}}) = ones(SMatrix{N,M,Float64})

#####################
## SMatrix methods ##
#####################

@pure size{S1,S2}(::Union{SMatrix{S1,S2},Type{SMatrix{S1,S2}}}) = (S1, S2)
@pure size{S1,S2,T}(::Type{SMatrix{S1,S2,T}}) = (S1, S2)
@pure size{S1,S2,T,L}(::Type{SMatrix{S1,S2,T,L}}) = (S1, S2)

function getindex(v::SMatrix, i::Integer)
    Base.@_inline_meta
    v.data[i]
end

@inline Tuple(v::SMatrix) = v.data

@inline function Base.unsafe_convert{N,M,T}(::Type{Ptr{T}}, m::SMatrix{N,M,T})
    Base.unsafe_convert(Ptr{T}, Base.data_pointer_from_objref(m))
end

macro SMatrix(ex)
    if !isa(ex, Expr)
        error("Bad input for @SMatrix")
    end

    if ex.head == :vect && length(ex.args) == 1 # 1 x 1
        return esc(Expr(:call, SMatrix{1, 1}, Expr(:tuple, ex.args[1])))
    elseif ex.head == :ref && length(ex.args) == 2 # typed, 1 x 1
        return esc(Expr(:call, Expr(:curly, :SMatrix, 1, 1, ex.args[1]), Expr(:tuple, ex.args[2])))
    elseif ex.head == :hcat # 1 x n
        s1 = 1
        s2 = length(ex.args)
        return esc(Expr(:call, SMatrix{s1, s2}, Expr(:tuple, ex.args...)))
    elseif ex.head == :typed_hcat # typed, 1 x n
        s1 = 1
        s2 = length(ex.args) - 1
        return esc(Expr(:call, Expr(:curly, :SMatrix, s1, s2, ex.args[1]), Expr(:tuple, ex.args[2:end]...)))
    elseif ex.head == :vcat
        if isa(ex.args[1], Expr) && ex.args[1].head == :row # n x m
            # Validate
            s1 = length(ex.args)
            s2s = map(i -> ((isa(ex.args[i], Expr) && ex.args[i].head == :row) ? length(ex.args[i].args) : 1), 1:s1)
            s2 = minimum(s2s)
            if maximum(s2s) != s2
                error("Rows must be of matching lengths")
            end

            exprs = [ex.args[i].args[j] for i = 1:s1, j = 1:s2]
            return esc(Expr(:call, SMatrix{s1, s2}, Expr(:tuple, exprs...)))
        else # n x 1
            return esc(Expr(:call, SMatrix{length(ex.args), 1}, Expr(:tuple, ex.args...)))
        end
    elseif ex.head == :typed_vcat
        if isa(ex.args[2], Expr) && ex.args[2].head == :row # typed, n x m
            # Validate
            s1 = length(ex.args) - 1
            s2s = map(i -> ((isa(ex.args[i+1], Expr) && ex.args[i+1].head == :row) ? length(ex.args[i+1].args) : 1), 1:s1)
            s2 = minimum(s2s)
            if maximum(s2s) != s2
                error("Rows must be of matching lengths")
            end

            exprs = [ex.args[i+1].args[j] for i = 1:s1, j = 1:s2]
            return esc(Expr(:call, Expr(:curly, :SMatrix,s1, s2, ex.args[1]), Expr(:tuple, exprs...)))
        else # typed, n x 1
            return esc(Expr(:call, Expr(:curly, :SMatrix, length(ex.args)-1, 1, ex.args[1]), Expr(:tuple, ex.args[2:end]...)))
        end
    elseif isa(ex, Expr) && ex.head == :comprehension
        if length(ex.args) != 1 || !isa(ex.args[1], Expr) || ex.args[1].head != :generator
            error("Expected generator in comprehension, e.g. [f(i,j) for i = 1:3, j = 1:3]")
        end
        ex = ex.args[1]
        if length(ex.args) != 3
            error("Use a 2-dimensional comprehension for @SMatrx")
        end

        rng1 = eval(current_module(), ex.args[2].args[2])
        rng2 = eval(current_module(), ex.args[3].args[2])
        f = gensym()
        f_expr = :($f = (($(ex.args[2].args[1]), $(ex.args[3].args[1])) -> $(ex.args[1])))
        exprs = [:($f($j1, $j2)) for j1 in rng1, j2 in rng2]

        return quote
            $(Expr(:meta, :inline))
            $(esc(f_expr))
            $(esc(Expr(:call, Expr(:curly, :SMatrix, length(rng1), length(rng2)), Expr(:tuple, exprs...))))
        end
    elseif isa(ex, Expr) && ex.head == :typed_comprehension
        if length(ex.args) != 2 || !isa(ex.args[2], Expr) || ex.args[2].head != :generator
            error("Expected generator in typed comprehension, e.g. Float64[f(i,j) for i = 1:3, j = 1:3]")
        end
        T = ex.args[1]
        ex = ex.args[2]
        if length(ex.args) != 3
            error("Use a 2-dimensional comprehension for @SMatrx")
        end

        rng1 = eval(current_module(), ex.args[2].args[2])
        rng2 = eval(current_module(), ex.args[3].args[2])
        f = gensym()
        f_expr = :($f = (($(ex.args[2].args[1]), $(ex.args[3].args[1])) -> $(ex.args[1])))
        exprs = [:($f($j1, $j2)) for j1 in rng1, j2 in rng2]

        return quote
            $(Expr(:meta, :inline))
            $(esc(f_expr))
            $(esc(Expr(:call, Expr(:curly, :SMatrix, length(rng1), length(rng2), T), Expr(:tuple, exprs...))))
        end
    elseif isa(ex, Expr) && ex.head == :call
        if ex.args[1] == :zeros || ex.args[1] == :ones || ex.args[1] == :rand || ex.args[1] == :randn
            if length(ex.args) == 3
                return quote
                    $(Expr(:meta, :inline))
                    $(ex.args[1])(SMatrix{$(esc(ex.args[2])),$(esc(ex.args[3]))})
                end
            elseif length(ex.args) == 4
                return quote
                    $(Expr(:meta, :inline))
                    $(ex.args[1])(SMatrix{$(esc(ex.args[3])), $(esc(ex.args[4])), $(esc(ex.args[2]))})
                end
            else
                error("@SMatrix expected a 2-dimensional array expression")
            end
        elseif ex.args[1] == :fill
            if length(ex.args) == 4
                return quote
                    $(Expr(:meta, :inline))
                    $(esc(ex.args[1]))($(esc(ex.args[2])), SMatrix{$(esc(ex.args[3])), $(esc(ex.args[4]))})
                end
            else
                error("@SMatrix expected a 2-dimensional array expression")
            end
        elseif ex.args[1] == :eye
            if length(ex.args) == 2
                return quote
                    $(Expr(:meta, :inline))
                    eye(SMatrix{$(esc(ex.args[2]))})
                end
            elseif length(ex.args) == 3
                # We need a branch, depending if the first argument is a type or a size.
                return quote
                    $(Expr(:meta, :inline))
                    if isa($(esc(ex.args[2])), DataType)
                        eye(SMatrix{$(esc(ex.args[3])), $(esc(ex.args[3])), $(esc(ex.args[2]))})
                    else
                        eye(SMatrix{$(esc(ex.args[2])), $(esc(ex.args[3]))})
                    end
                end
            elseif length(ex.args) == 4
                return quote
                    $(Expr(:meta, :inline))
                    eye(SMatrix{$(esc(ex.args[3])), $(esc(ex.args[4])), $(esc(ex.args[2]))})
                end
            else
                error("Bad eye() expression for @SMatrix")
            end
        else
            error("@SMatrix only supports the zeros(), ones(), fill(), rand(), randn() and eye() functions.")
        end
    else
        error("Bad input for @SMatrix")
    end
end


text:
# ------------------------------------------------------------------------------------------------ #
#          Run a statistical test to compare the distributions of features in two matrices         #
# ------------------------------------------------------------------------------------------------ #
function statfeaturefilter(F₁, F₂; tail=:firstlarger, test=TailedTestOfVariance, alpha=0.05)

    p = testFeatureDistributions(F₁, F₂, test; tail=tail)
    issig = p .< alpha
    strs = deepcopy(p)
    strs[issig] = 10.0.^(ceil.(log10.(p[issig])))
    strs[.!issig] = round.(strs[.!issig], sigdigits=2)
    strs = string.(strs)
    strs[issig] = "<".*strs[issig]
    ann = [text("$(strs[i])", [:crimson, :black][issig[i]+1], 8, rotation=90, valign=:bottom) for i ∈ 1:length(p)] # :forestgreen
    plt = featureviolin(F₁, F₂, rightannotations=ann)
end


text:
using GoogleCloud

CREDENTIALS = nothing

function init(filepath)
    global CREDENTIALS
    CREDENTIALS = JSONCredentials(filepath)
end

function projectid()
    global CREDENTIALS
    isnothing(CREDENTIALS) && throw("Run JuliaFirestore.init(filepath) to load credentials first")
    CREDENTIALS.project_id
end

function get_token()
    global CREDENTIALS
    isnothing(CREDENTIALS) && throw("Run JuliaFirestore.init(filepath) to load credentials first")
    session = GoogleSession(CREDENTIALS, ["datastore"])
    GoogleCloud.authorize(session)
end

function auth_header()
    token = get_token()
    "Authorization" => "$(token[:token_type]) $(token[:access_token])"
end


text:
using SortMark

df = make_df()

skip = (df.Type .<: Union{UInt64,UInt128}) .&
    (df.source_key .== :small_positive)
compute!(df[(!).(skip), :])
stat!(df)
SortMark

max_stars_repo_name: LilithHafner/SortMark

max_issues_repo_name: LilithHafner/SortMark


text:
# -------------------------------------------------------------------------------------------------------------------------------------------------------------
# Script producing raw data for BE
# -------------------------------------------------------------------------------------------------------------------------------------------------------------

using DataFrames, CSV, XLSX, LinearAlgebra, Statistics, Dates

# other scripts
dir_home = "x:/VIVES/1-Personal/Florian/git/2022_01_comext/src/"
include(dir_home * "functions.jl")

# location of data input/output (io)
dir_io = "C:/Users/u0148308/data/comext/"
dir_dropbox = "C:/Users/u0148308/Dropbox/BREXIT/"

# -------------------------------------------------------------------------------------------------------------------------------------------------------------
# Download raw data from Comext Bulk Download Facility
# -------------------------------------------------------------------------------------------------------------------------------------------------------------

url = "https://ec.europa.eu/eurostat/estat-navtree-portlet-prod/BulkDownloadListing?sort=1&file=comext%2FCOMEXT_DATA%2FPRODUCTS%2F"

years = string.(2001:2021)
months = lpad.(1:12, 2, '0')

for i in years
    for j in months

        id = "full" * i * j * ".7z"

        if !isfile(dir_dropbox * "rawdata/comext/zipped/" * id)
            download(url * id, dir_dropbox * "rawdata/comext/zipped/" * id)
        else
            println(" ✓ The zipped file for $j/$i has already been downloaded.")
        end

    end
end

# -------------------------------------------------------------------------------------------------------------------------------------------------------------
# Unzip raw data
# -------------------------------------------------------------------------------------------------------------------------------------------------------------

# trying to automate un-zipping, unfortunately failed => NEED TO MANUALLY UNZIPP to folder dir_io!
# run(`cmd /c set PATH=%PATH% ';' "C:\\Program Files\\7-Zip\\" echo %PATH% 7z`)
# run(`cmd /c cd "C:\\Users\\u0148308\\data\\raw\\"`)
# run(`cmd /c 7z e a.7z`)

# manually

# -------------------------------------------------------------------------------------------------------------------------------------------------------------
# Initial cleaning, subset and export data for BE as DECLARANT
# -------------------------------------------------------------------------------------------------------------------------------------------------------------

# EU27
ctrys_EU27 = ["AT", "BE", "BG", "CY", "CZ", "DE", "DK", "EE", "ES", "FI", "FR", "GR", "HR", "HU",
             "IE", "IT", "LT", "LU", "LV", "MT", "NL", "PL", "PT", "RO", "SE", "SI", "SK"]

for i in years
    for j in months

        path = dir_dropbox * "rawdata/comext/BE/" * "full" * i * j * "_BE" * ".csv"

        if isfile(path)
            println(" ✓ Data for BE in $j/$i already exists! \n")
        else
            # import and clean data
            df = initial_cleaning(i, j)

            # subset to include only BE as DECLARANT
            subset!(df, :DECLARANT_ISO => ByRow(x -> x == "BE"))

            # create EU/WORLD as PARTNER
            df = append_WORLD(df)
            df = append_EU(df, ctrys_EU27)

            # export to Dropbox
            CSV.write(path, df)

            println(" ✓ Data for BE in $j/$i has been successfully exported to Dropbox. \n")
        end

    end
end

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


text:
module MosaicViews

using PaddedViews
using OffsetArrays

export

    MosaicView,
    mosaicview

"""
    MosaicView(A::AbstractArray)

Create a two dimensional "view" of the three or four dimensional
array `A`. The resulting `MosaicView` will display the data in
`A` such that it emulates using `vcat` for all elements in the
third dimension of `A`, and `hcat` for all elements in the fourth
dimension of `A`.

For example, if `size(A)` is `(2,3,4)`, then the resulting
`MosaicView` will have the size `(2*4,3)` which is `(8,3)`.
Alternatively, if `size(A)` is `(2,3,4,5)`, then the resulting
size will be `(2*4,3*5)` which is `(8,15)`.

Another way to think about this is that `MosaicView` creates a
mosaic of all the individual matrices enumerated in the third
(and optionally fourth) dimension of the given 3D or 4D array
`A`. This can be especially useful for creating a single
composite image from a set of equally sized images.

```@jldoctest
julia> using MosaicViews

julia> A = [(k+1)*l-1 for i in 1:2, j in 1:3, k in 1:2, l in 1:2]
2×3×2×2 Array{Int64,4}:
[:, :, 1, 1] =
 1  1  1
 1  1  1

[:, :, 2, 1] =
 2  2  2
 2  2  2

[:, :, 1, 2] =
 3  3  3
 3  3  3

[:, :, 2, 2] =
 5  5  5
 5  5  5

julia> MosaicView(A)
4×6 MosaicViews.MosaicView{Int64,4,Array{Int64,4}}:
 1  1  1  3  3  3
 1  1  1  3  3  3
 2  2  2  5  5  5
 2  2  2  5  5  5
```
"""
struct MosaicView{T,N,A<:AbstractArray{T,N}} <: AbstractArray{T,2}
    parent::A
    dims::Tuple{Int,Int}
    pdims::NTuple{N,Int}

    function MosaicView{T,N}(A::AbstractArray{T,N}, dims) where {T,N}
        1 <= N <= 4 || throw(ArgumentError("The given array must have dimensionality N=3 or N=4"))
        # unless we store the axes in the struct, we can't support offset indices
        # N=2 is a special case that we can provide a specialization on `axes`
        # but for consistency with cases of other dimensions, disable it as well
        require_one_based_indexing(A)
        new{T,N,typeof(A)}(A, dims, size(A))
    end
end

function MosaicView(A::AbstractArray{T,N}) where {T,N}
    # vectors are lifted to 2D matrix, 3d/4d arrays are reshaped to 2D matrix
    dims = (size(A,1) * size(A,3), size(A,2) * size(A,4))
    MosaicView{T,N}(A, dims)
end

Base.parent(mv::MosaicView) = mv.parent
Base.size(mv::MosaicView) = mv.dims

# fallback for 1d/2d case
@inline Base.getindex(mv::MosaicView, ind::Int...) = mv.parent[ind...]

@inline function Base.getindex(mv::MosaicView{T,3,A}, i::Int, j::Int) where {T,A}
    @boundscheck checkbounds(mv, i, j)
    pdims = mv.pdims
    parent = mv.parent
    idx1 = (i-1) % pdims[1] + 1
    idx2 = (j-1) % pdims[2] + 1
    idx3 = (i-1) ÷ pdims[1] + 1
    @inbounds res = parent[idx1, idx2, idx3]
    res
end

@inline function Base.getindex(mv::MosaicView{T,4,A}, i::Int, j::Int) where {T,A}
    @boundscheck checkbounds(mv, i, j)
    pdims = mv.pdims
    parent = mv.parent
    idx1 = (i-1) % pdims[1] + 1
    idx2 = (j-1) % pdims[2] + 1
    idx3 = (i-1) ÷ pdims[1] + 1
    idx4 = (j-1) ÷ pdims[2] + 1
    @inbounds res = parent[idx1, idx2, idx3, idx4]
    res
end

"""
    mosaicview(A::AbstractArray;
               [fillvalue=<zero unit>], [npad=0],
               [nrow], [ncol], [rowmajor=false],
               [center=true]) -> MosaicView
    mosaicview(As::AbstractArray...; kwargs...)
    mosaicview(As::Union{Tuple, AbstractVector}; kwargs...)

Create a two dimensional "view" from array `A` or a list of arrays `As`.

The resulting [`MosaicView`](@ref) will display all the matrix
slices of the first two dimensions of `A` arranged as a single
large mosaic (in the form of a matrix).

If multiple arrays in passed, they'll will be center-padded to a common
size, and then be concatenated to create a higher dimensional array.

# Arguments

In contrast to using the constructor of [`MosaicView`](@ref)
directly, the function `mosaicview` also allows for a couple of
convenience keywords. A typical use case would be to create an
image mosaic from a set of input images.

- The parameter `fillvalue` defines the value that
  that should be used for empty space. This can be padding caused
  by `npad`, or empty mosaic tiles in case the number of matrix
  slices in `A` is smaller than `nrow*ncol`.

- The parameter `npad` defines the empty padding space between
  adjacent mosaic tiles. This can be especially useful if the
  individual tiles (i.e. matrix slices in `A`) are images that
  should be visually separated by some grid lines.

- The parameters `nrow` and `ncol` can be used to choose the
  number of tiles in row and/or column direction the mosaic should
  be arranged in. Note that it suffices to specify one of the
  two parameters, as the other one can be inferred accordingly.
  The default in case none of the two are specified is `nrow = size(A,3)`.

- If `rowmajor` is set to `true`, then the slices will be
  arranged left-to-right-top-to-bottom, instead of
  top-to-bottom-left-to-right (default). The layout only differs
  in non-trivial cases, i.e., when `nrow != 1` and `ncol != 1`.

- If `center` is set to `true`, then the padded arrays will be shifted
  to the center instead of in the top-left corner (default). This
  parameter is only useful when arrays are of different sizes.

!!! tip
    This function is not type stable and should only be used if
    performance is not a priority. To achieve optimized performance,
    you need to manually construct a [`MosaicView`](@ref).

# Examples

The simplest usage is to `cat` two arrays of the same dimension.

```julia-repl
julia> A1 = fill(1, 3, 1)
3×1 Array{Int64,2}:
 1
 1
 1

julia> A2 = fill(2, 1, 3)
1×3 Array{Int64,2}:
 2  2  2

julia> mosaicview(A1, A2)
6×3 MosaicView{Int64,4, ...}:
 0  1  0
 0  1  0
 0  1  0
 0  0  0
 2  2  2
 0  0  0

julia> mosaicview(A1, A2; center=false)
 6×3 MosaicView{Int64,4, ...}:
  1  0  0
  1  0  0
  1  0  0
  2  2  2
  0  0  0
  0  0  0
```

Other keyword arguments can be useful to get a nice looking results.

```julia-repl
julia> using MosaicViews

julia> A = [k for i in 1:2, j in 1:3, k in 1:5]
2×3×5 Array{Int64,3}:
[:, :, 1] =
 1  1  1
 1  1  1

[:, :, 2] =
 2  2  2
 2  2  2

[:, :, 3] =
 3  3  3
 3  3  3

[:, :, 4] =
 4  4  4
 4  4  4

[:, :, 5] =
 5  5  5
 5  5  5

julia> mosaicview(A, ncol=2)
6×6 MosaicViews.MosaicView{Int64,4,...}:
 1  1  1  4  4  4
 1  1  1  4  4  4
 2  2  2  5  5  5
 2  2  2  5  5  5
 3  3  3  0  0  0
 3  3  3  0  0  0

julia> mosaicview(A, nrow=2)
4×9 MosaicViews.MosaicView{Int64,4,...}:
 1  1  1  3  3  3  5  5  5
 1  1  1  3  3  3  5  5  5
 2  2  2  4  4  4  0  0  0
 2  2  2  4  4  4  0  0  0

julia> mosaicview(A, nrow=2, rowmajor=true)
4×9 MosaicViews.MosaicView{Int64,4,...}:
 1  1  1  2  2  2  3  3  3
 1  1  1  2  2  2  3  3  3
 4  4  4  5  5  5  0  0  0
 4  4  4  5  5  5  0  0  0

julia> mosaicview(A, nrow=2, npad=1, rowmajor=true)
5×11 MosaicViews.MosaicView{Int64,4,...}:
 1  1  1  0  2  2  2  0  3  3  3
 1  1  1  0  2  2  2  0  3  3  3
 0  0  0  0  0  0  0  0  0  0  0
 4  4  4  0  5  5  5  0  0  0  0
 4  4  4  0  5  5  5  0  0  0  0

julia> mosaicview(A, fillvalue=-1, nrow=2, npad=1, rowmajor=true)
5×11 MosaicViews.MosaicView{Int64,4,...}:
  1   1   1  -1   2   2   2  -1   3   3   3
  1   1   1  -1   2   2   2  -1   3   3   3
 -1  -1  -1  -1  -1  -1  -1  -1  -1  -1  -1
  4   4   4  -1   5   5   5  -1  -1  -1  -1
  4   4   4  -1   5   5   5  -1  -1  -1  -1
```
"""
function mosaicview(A::AbstractArray{T,3};
                    fillvalue = zero(T),
                    npad = 0,
                    nrow = -1,
                    ncol = -1,
                    rowmajor = false,
                    kwargs...) where T
    nrow == -1 || nrow > 0 || throw(ArgumentError("The parameter \"nrow\" must be greater than 0"))
    ncol == -1 || ncol > 0 || throw(ArgumentError("The parameter \"ncol\" must be greater than 0"))
    npad >= 0 || throw(ArgumentError("The parameter \"npad\" must be greater than or equal to 0"))
    ntile = size(A,3)
    ntile_ceil = ntile # ntile need not be integer divideable
    if nrow == -1 && ncol == -1
        # automatically choose nrow to reflect what MosaicView does
        nrow = ntile
        ncol = 1
    elseif nrow == -1
        # compute nrow based on ncol
        nrow = ceil(Int, ntile / ncol)
        ntile_ceil = nrow * ncol
    elseif ncol == -1
        # compute ncol based on nrow
        ncol = ceil(Int, ntile / nrow)
        ntile_ceil = nrow * ncol
    else
        # accept nrow and ncol as is if it covers at least all
        # existing tiles
        ntile_ceil = nrow * ncol
        ntile_ceil < ntile && throw(ArgumentError("The product of the parameters \"ncol\" (value: $ncol) and \"nrow\" (value: $nrow) must be greater than or equal to $ntile"))
    end
    # we pad size(A,3) to nrow*ncol. we also pad the first two
    # dimensions according to npad. think of this as border
    # between tiles (useful for images)
    pad_dims = (size(A,1) + npad, size(A,2) + npad, ntile_ceil)
    A_pad = PaddedView(T(fillvalue), A, pad_dims)
    # next we reshape the image such that it reflects the
    # specified nrow and ncol
    A_new = if !rowmajor
        res_dims = (size(A_pad,1), size(A_pad,2), nrow, ncol)
        reshape(A_pad, res_dims)
    else
        # same as above but we additionally permute dimensions
        # to mimic row first layout for the tiles. this is useful
        # for images since user often reads these from left-to-right
        # before top-to-bottom. (note the swap of "ncol" and "nrow")
        res_dims = (size(A_pad,1), size(A_pad,2), ncol, nrow)
        A_tp = reshape(A_pad, res_dims)
        PermutedDimsArray(A_tp, (1, 2, 4, 3))
    end
    # decrease size of the resulting MosaicView by npad to not have
    # a border on the right side and bottom side of the final mosaic.
    dims = (size(A_new,1) * size(A_new,3) - npad, size(A_new,2) * size(A_new,4) - npad)
    MosaicView{T,4}(A_new, dims)
end

function mosaicview(A::AbstractArray{T,N};
                    nrow = -1,
                    ncol = -1,
                    kwargs...) where {T,N}
    # if neither nrow nor ncol is provided then automatically choose
    # nrow and ncol to reflect what MosaicView does (i.e. use size)
    if nrow == -1 && ncol == -1
        nrow = size(A, 3)
        # ncol = size(A, 4)
    end
    mosaicview(reshape(A, (size(A,1), size(A,2), :));
               nrow=nrow, ncol=ncol, kwargs...)
end

mosaicview(As::AbstractArray...; kwargs...) = mosaicview(As; kwargs...)

function mosaicview(As::AbstractVector{T};
                    fillvalue=zero(eltype(first(As))),
                    center=true,
                    kwargs...) where {T <: AbstractArray}
    length(As) == 0 && throw(ArgumentError("The given vector should not be empty"))
    length(unique(ndims.(As))) != 1 && throw(ArgumentError("All arrays should have the same dimension"))
    N = ndims(first(As))
    mosaicview(_padded_cat(As; center=center, fillvalue=fillvalue, dims=max(3, N+1));
               fillvalue=fillvalue, kwargs...)
end

function mosaicview(As::Tuple;
                    fillvalue=zero(eltype(first(As))),
                    center=true,
                    kwargs...)
    length(As) == 0 && throw(ArgumentError("The given tuple should not be empty"))
    length(unique(ndims.(As))) != 1 && throw(ArgumentError("All arrays should have the same dimension"))
    N = ndims(first(As))
    mosaicview(_padded_cat(As; center=center, fillvalue=fillvalue, dims=max(3, N+1));
               fillvalue=fillvalue, kwargs...)
end

function _padded_cat(imgs; center, fillvalue, dims)
    # reduce(cat, imgs) would indeed make the whole pipeline more eagerly
    # and thus allocates more memory
    # TODO: inefficient when there're too many images, e.g., 512
    if length(imgs) > 300
        msg = "It's quite slow to visualize a tuple/list of $(length(imgs)) images"
        msg *= "\nyou might want to manually `cat` them into one large array first."
        @warn msg
    end
    has_offsets = any(Base.has_offset_axes.(imgs))
    if !has_offsets && length(unique(map(axes, imgs))) == 1
        return cat(imgs...; dims=dims)
    else
        if !has_offsets && !center
            # in this case the index of PaddedView starts from 1
            # hence we can directly pass them into `cat`
            cat(paddedviews(fillvalue, imgs...)...; dims=dims)
        else
            # TODO: it's unidentified but this version allocates more memory
            #       and become slower (~1.5X) than a direct `cat`
            pad_fn = center ? sym_paddedviews : paddedviews
            reduce(pad_fn(fillvalue, imgs...)) do x, y
                x = OffsetArray(x, 1 .- first.(axes(x)))
                y = OffsetArray(y, 1 .- first.(axes(y)))
                cat(x, y; dims=dims)
            end
        end
    end
end

### compat
if VERSION <= v"1.1"
    require_one_based_indexing(A...) = !Base.has_offset_axes(A...) || throw(ArgumentError("offset arrays are not supported but got an array with index other than 1"))
else
    require_one_based_indexing = Base.require_one_based_indexing
end
### deprecations

@deprecate mosaicview(A::AbstractArray, fillvalue; kwargs...) mosaicview(A; fillvalue=fillvalue, kwargs...)

end # module


text:
struct Adams_Bashforth3rdO{Adaptative} <: AbstractScalarTimeStep{Adaptative,2}
    fm1::ScalarField{Float64,3,2,false} #Store latest step
    fm2::ScalarField{Float64,3,2,false} #Store 2 steps before
    At::Base.RefValue{Float64} # 23*dt/12 for constant time stepping
    Bt::Base.RefValue{Float64} # -16*dt/12 for constant time stepping
    Ct::Base.RefValue{Float64} # 5*dt/12 for constant time stepping
    dt::Base.RefValue{Float64}
    dt2::Base.RefValue{Float64}
    dt3::Base.RefValue{Float64}
    iteration::Base.RefValue{Int}
end

function Adams_Bashforth3rdO(adp,indt)
    at = Ref(23*indt/12)
    bt = Ref(-16*indt/12)
    ct = Ref(5*indt/12)
    dt = Ref(indt)
    dt2 = Ref(indt)
    dt3 = Ref(indt)
    fm1 = ScalarField{Float64}((NRX,NY,NZ),(LX,LY,LZ))
    fm2 = ScalarField{Float64}((NRX,NY,NZ),(LX,LY,LZ))
    return Adams_Bashforth3rdO{adp}(fm1,fm2,at,bt,ct,dt,dt2,dt3,Ref(0))
end

get_dt(t::Adams_Bashforth3rdO) =
    getindex(t.dt)

function initialize!(t::Adams_Bashforth3rdO,rhs::AbstractArray,vis,s::AbstractSimulation)
    mycopy!(data(t.fm1), rhs)
    mycopy!(data(t.fm2), data(t.fm1))
    setindex!(t.dt,get_dt(s))
    setindex!(t.dt2,t.dt[])
    setindex!(t.dt3,t.dt[])
    return nothing
end

function set_dt!(t::Adams_Bashforth3rdO{true},dt::Real)
    if t.iteration[] == 0
        setindex!(t.dt3,dt/2)
        setindex!(t.dt2,dt/2)
        setindex!(t.dt,dt/2)
    else
        setindex!(t.dt3,t.dt2[])
        setindex!(t.dt2,t.dt[])
        setindex!(t.dt,dt)
    end
end

function get_At(t::Adams_Bashforth3rdO{true})
    dt = t.dt[]
    dt2 = t.dt2[]
    dt3 = t.dt3[]

    At = 1. + (dt*(2*dt + 6*dt2 + 3*dt3))/(6*dt2*(dt2 + dt3))
    At *= dt
    return At
end

get_At(t::Type{Adams_Bashforth3rdO{false}}) =
    23*get_dt(t)/12

@inline get_At(t::A) where {A<:Adams_Bashforth3rdO{false}} =
    get_At(A)

function get_Bt(t::Adams_Bashforth3rdO{true})
    dt = t.dt[]
    dt2 = t.dt2[]
    dt3 = t.dt3[]

    Bt = -dt*dt*(2*dt + 3*dt2 + 3*dt3)/(6*dt2*dt3*(dt+dt2))
    Bt *= (dt+dt2)
    return Bt
end

get_Bt(t::Type{Adams_Bashforth3rdO{false}}) =
    -16*get_dt(t)/12

@inline get_Bt(t::A) where {A<:Adams_Bashforth3rdO{false}} =
    get_Bt(A)

function get_Ct(t::Adams_Bashforth3rdO{true})
    dt = t.dt[]
    dt2 = t.dt2[]
    dt3 = t.dt3[]

    Ct = dt*dt*(2*dt + 3*dt2)/(6*dt3*(dt + dt2 + dt3)*(dt2 + dt3))
    Ct*=(dt+dt2+dt3)
    return Ct
end

get_Ct(t::Type{Adams_Bashforth3rdO{false}}) =
    5*get_det(t)/12

@inline get_Ct(t::A) where {A<:Adams_Bashforth3rdO{false}} =
    get_Ct(A)

function set_coefficients!(t::Adams_Bashforth3rdO,rhs)
    i = t.iteration[] += 1
    i == 2 && fix_fm2!(t.fm2,t.fm1,rhs,t.dt2[],t.dt3[])
end

@par function (f::Adams_Bashforth3rdO)(ρ::AbstractArray{<:Complex,3},ρrhs::AbstractArray{<:Complex,3}, s::@par(AbstractSimulation))
    set_coefficients!(f,ρrhs)
    @mthreads for ind in ZYRANGE
        jj,kk = Tuple(ind)
        k = ZRANGE[kk]
        j = YRANGE[jj]
        _tAdams_Bashforth3rdO!(j,k, ρ,ρrhs,f.fm1,f.fm2,f,s)
    end
    return nothing
end

@inline @par function _tAdams_Bashforth3rdO!(j::Integer,k::Integer, u::AbstractArray{Complex{Float64},3}, rhs::AbstractArray, rm1::AbstractArray, rm2::AbstractArray, f, s::@par(AbstractSimulation))
@inbounds begin
    At = get_At(f)
    Bt = get_Bt(f)
    Ct = get_Ct(f)
    @inbounds begin
        @msimd for i in XRANGE
            #u[i] += dt12*(23*rhs[i] - 16rm1[i] + 5rm2[i])
            #u[i,j,k] = muladd(muladd(23, rhs[i,j,k], muladd(-16, rm1[i,jj,kk], 5rm2[i,jj,kk])), dt12, u[i,j,k])
            u[i,j,k] = muladd(At, rhs[i,j,k], muladd(Bt, rm1[i,j,k], muladd(Ct, rm2[i,j,k], u[i,j,k])))
            rm2[i,j,k] = rm1[i,j,k]
            rm1[i,j,k] = rhs[i,j,k]
        end
    end
end
end


text:
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#
# Description
# ==============================================================================
#
#   Constants use to compute the moon position.
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#
# References
# ==============================================================================
#
#   [1] Meeus, J (1998). Astronomical algorithms. Willmann-Bell, Inc, Richmond,
#       VA.
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #

# Those tables were obtained from [1, p. 339]. They are used to compute the Moon
# position vector.
const _tab_47a  = [
    0  0  1  0 6288774 -20905355
    2  0 -1  0 1274027  -3699111
    2  0  0  0  658314  -2955968
    0  0  2  0  213618   -569925
    0  1  0  0 -185116     48888
    0  0  0  2 -114332     -3149
    2  0 -2  0   58793    246158
    2 -1 -1  0   57066   -152138
    2  0  1  0   53322   -170733
    2 -1  0  0   45758   -204586
    0  1 -1  0  -40923   -129620
    1  0  0  0  -34720    108743
    0  1  1  0  -30383    104755
    2  0  0 -2   15327     10321
    0  0  1  2  -12528         0
    0  0  1 -2   10980     79661
    4  0 -1  0   10675    -34782
    0  0  3  0   10034    -23210
    4  0 -2  0    8548    -21636
    2  1 -1  0   -7888     24208
    2  1  0  0   -6766     30824
    1  0 -1  0   -5163     -8379
    1  1  0  0    4987    -16675
    2 -1  1  0    4036    -12831
    2  0  2  0    3994    -10445
    4  0  0  0    3861    -11650
    2  0 -3  0    3665     14403
    0  1 -2  0   -2689     -7003
    2  0 -1  2   -2602         0
    2 -1 -2  0    2390     10056
    1  0  1  0   -2348      6322
    2 -2  0  0    2236     -9884
    0  1  2  0   -2120      5751
    0  2  0  0   -2069         0
    2 -2 -1  0    2048     -4950
    2  0  1 -2   -1773      4130
    2  0  0  2   -1595         0
    4 -1 -1  0    1215     -3958
    0  0  2  2   -1110         0
    3  0 -1  0    -892      3258
    2  1  1  0    -810      2616
    4 -1 -2  0     759     -1897
    0  2 -1  0    -713     -2117
    2  2 -1  0    -700      2354
    2  1 -2  0     691         0
    2 -1  0 -2     596         0
    4  0  1  0     549     -1423
    0  0  4  0     537     -1117
    4 -1  0  0     520     -1571
    1  0 -2  0    -487     -1739
    2  1  0 -2    -399         0
    0  0  2 -2    -381     -4421
    1  1  1  0     351         0
    3  0 -2  0    -340         0
    4  0 -3  0     330         0
    2 -1  2  0     327         0
    0  2  1  0    -323      1165
    1  1 -1  0     299         0
    2  0  3  0     294         0
    2  0 -1 -2       0     8752
]'

const _tab_47b = [
    0  0  0  1 5128122
    0  0  1  1  280602
    0  0  1 -1  277693
    2  0  0 -1  173237
    2  0 -1  1   55413
    2  0 -1 -1   46271
    2  0  0  1   32573
    0  0  2  1   17198
    2  0  1 -1    9266
    0  0  2 -1    8822
    2 -1  0 -1    8216
    2  0 -2 -1    4324
    2  0  1  1    4200
    2  1  0 -1   -3359
    2 -1 -1  1    2463
    2 -1  0  1    2211
    2 -1 -1 -1    2065
    0  1 -1 -1   -1870
    4  0 -1 -1    1828
    0  1  0  1   -1794
    0  0  0  3   -1749
    0  1 -1  1   -1565
    1  0  0  1   -1491
    0  1  1  1   -1475
    0  1  1 -1   -1410
    0  1  0 -1   -1344
    1  0  0 -1   -1335
    0  0  3  1    1107
    4  0  0 -1    1021
    4  0 -1  1     833
    0  0  1 -3     777
    4  0 -2  1     671
    2  0  0 -3     607
    2  0  2 -1     596
    2 -1  1 -1     491
    2  0 -2  1    -451
    0  0  3 -1     439
    2  0  2  1     422
    2  0 -3 -1     421
    2  1 -1  1    -366
    2  1  0  1    -351
    4  0  0  1     331
    2 -1  1  1     315
    2 -2  0 -1     302
    0  0  1  3    -283
    2  1  1 -1    -229
    1  1  0 -1     223
    1  1  0  1     223
    0  1 -2 -1    -220
    2  1 -1 -1    -220
    1  0  1  1    -185
    2 -1 -2 -1     181
    0  1  2  1    -177
    4  0 -2 -1     176
    4 -1 -1 -1     166
    1  0  1 -1    -164
    4  0  1 -1     132
    1  0 -1 -1    -119
    4 -1  0 -1     115
    2 -2  0  1     107
]'


text:
struct CommissionReport
  execId::String
  commission::Float64
  currency::String
  realizedPNL::Union{Float64,Nothing}
  yield::Union{Float64,Nothing}
  yieldRedemptionDate::Union{Int,Nothing}     # yyyymmdd format TODO: -> String or Date
end

struct ContractDescription
  contract::Contract
  derivativeSecTypes::Vector{String}
end

mutable struct ContractDetails
  contract::Contract
  marketName::String
  minTick::Float64
  orderTypes::String
  validExchanges::String
  priceMagnifier::Int
  underConId::Int
  longName::String
  contractMonth::String
  industry::String
  category::String
  subcategory::String
  timeZoneId::String
  tradingHours::String
  liquidHours::String
  evRule::String
  evMultiplier::Union{Float64,Nothing}
  mdSizeMultiplier::Union{Int,Nothing}
  aggGroup::Union{Int,Nothing}
  underSymbol::String
  underSecType::String
  marketRuleIds::String
  realExpirationDate::String
  lastTradeTime::String
  stockType::String
  secIdList::NamedTuple
  cusip::String
  ratings::String
  descAppend::String
  bondType::String
  couponType::String
  callable::Bool
  putable::Bool
  coupon::Float64
  convertible::Bool
  maturity::String
  issueDate::String
  nextOptionDate::String
  nextOptionType::String
  nextOptionPartial::Bool
  notes::String
end
ContractDetails() = ContractDetails(Contract(), ns, 0.0, ns, ns, 0, 0, fill(ns, 9)..., 0,
                                    nothing, nothing, fill(ns, 6)..., NamedTuple(),
                                    fill(ns, 5)..., false, false, 0, false, ns, ns, ns, ns,
                                    false, ns)

struct Execution
  orderId::Union{Int,Nothing}     # Moved up from position 10 for convenience when constructing
  execId::String
  time::String
  acctNumber::String
  exchange::String
  side::String
  shares::Float64
  price::Float64
  permId::Int
  clientId::Int
  liquidation::Int
  cumQty::Float64
  avgPrice::Float64
  orderRef::String
  evRule::String
  evMultiplier::Union{Float64,Nothing}
  modelCode::String
  lastLiquidity::Int
end

struct FamilyCode
  accountID::String
  familyCodeStr::String
end

struct OrderState
  status::String
  initMarginBefore::String
  maintMarginBefore::String
  equityWithLoanBefore::String
  initMarginChange::String
  maintMarginChange::String
  equityWithLoanChange::String
  initMarginAfter::String
  maintMarginAfter::String
  equityWithLoanAfter::String
  commission::Union{Float64,Nothing}
  minCommission::Union{Float64,Nothing}
  maxCommission::Union{Float64,Nothing}
  commissionCurrency::String
  warningText::String
  completedTime::String
  completedStatus::String
end

TickAttrib =       NamedTuple{(:canAutoExecute, :pastLimit, :preOpen),NTuple{3,Bool}}
TickAttribLast =   NamedTuple{(:pastLimit, :unreported),NTuple{2,Bool}}
TickAttribBidAsk = NamedTuple{(:bidPastLow, :askPastHigh),NTuple{2,Bool}}


text:
using JuMP, EAGO

                     m = Model()

                     EAGO.register_eago_operators!(m)

                     @variable(m, -1 <= x[i=1:4] <= 1)

                     @variable(m, -2.9922365365158297 <= q <= 1.3959867501971397)

                     add_NL_constraint(m, :(softplus(-0.7172866311917714 + 0.8530094465195428*softplus(0.49975691401319544 + 0.817850285997455*softplus(-0.9044530867886094 + -0.3504940249833268*softplus(0.9239332442221326 + 0.41157080299867177*$(x[1]) + 0.6856618842893334*$(x[2]) + 0.23250449051531152*$(x[3]) + -0.8171094215290706*$(x[4])) + -0.925387475816613*softplus(-0.1441575648301785 + -0.7324773919132506*$(x[1]) + 0.1666248349273194*$(x[2]) + 0.5034013066229974*$(x[3]) + -0.3351866390106042*$(x[4]))) + -0.15011340339301915*softplus(-0.49866226232151334 + 0.31623027699891626*softplus(0.9239332442221326 + 0.41157080299867177*$(x[1]) + 0.6856618842893334*$(x[2]) + 0.23250449051531152*$(x[3]) + -0.8171094215290706*$(x[4])) + 0.1823838704363938*softplus(-0.1441575648301785 + -0.7324773919132506*$(x[1]) + 0.1666248349273194*$(x[2]) + 0.5034013066229974*$(x[3]) + -0.3351866390106042*$(x[4])))) + -0.2302663292873568*softplus(0.6056562038214484 + 0.017530899400252764*softplus(-0.9044530867886094 + -0.3504940249833268*softplus(0.9239332442221326 + 0.41157080299867177*$(x[1]) + 0.6856618842893334*$(x[2]) + 0.23250449051531152*$(x[3]) + -0.8171094215290706*$(x[4])) + -0.925387475816613*softplus(-0.1441575648301785 + -0.7324773919132506*$(x[1]) + 0.1666248349273194*$(x[2]) + 0.5034013066229974*$(x[3]) + -0.3351866390106042*$(x[4]))) + 0.9141090581590978*softplus(-0.49866226232151334 + 0.31623027699891626*softplus(0.9239332442221326 + 0.41157080299867177*$(x[1]) + 0.6856618842893334*$(x[2]) + 0.23250449051531152*$(x[3]) + -0.8171094215290706*$(x[4])) + 0.1823838704363938*softplus(-0.1441575648301785 + -0.7324773919132506*$(x[1]) + 0.1666248349273194*$(x[2]) + 0.5034013066229974*$(x[3]) + -0.3351866390106042*$(x[4]))))) + softplus(0.11802092343800519 + -0.23198020992542823*softplus(0.49975691401319544 + 0.817850285997455*softplus(-0.9044530867886094 + -0.3504940249833268*softplus(0.9239332442221326 + 0.41157080299867177*$(x[1]) + 0.6856618842893334*$(x[2]) + 0.23250449051531152*$(x[3]) + -0.8171094215290706*$(x[4])) + -0.925387475816613*softplus(-0.1441575648301785 + -0.7324773919132506*$(x[1]) + 0.1666248349273194*$(x[2]) + 0.5034013066229974*$(x[3]) + -0.3351866390106042*$(x[4]))) + -0.15011340339301915*softplus(-0.49866226232151334 + 0.31623027699891626*softplus(0.9239332442221326 + 0.41157080299867177*$(x[1]) + 0.6856618842893334*$(x[2]) + 0.23250449051531152*$(x[3]) + -0.8171094215290706*$(x[4])) + 0.1823838704363938*softplus(-0.1441575648301785 + -0.7324773919132506*$(x[1]) + 0.1666248349273194*$(x[2]) + 0.5034013066229974*$(x[3]) + -0.3351866390106042*$(x[4])))) + 0.3694431087223533*softplus(0.6056562038214484 + 0.017530899400252764*softplus(-0.9044530867886094 + -0.3504940249833268*softplus(0.9239332442221326 + 0.41157080299867177*$(x[1]) + 0.6856618842893334*$(x[2]) + 0.23250449051531152*$(x[3]) + -0.8171094215290706*$(x[4])) + -0.925387475816613*softplus(-0.1441575648301785 + -0.7324773919132506*$(x[1]) + 0.1666248349273194*$(x[2]) + 0.5034013066229974*$(x[3]) + -0.3351866390106042*$(x[4]))) + 0.9141090581590978*softplus(-0.49866226232151334 + 0.31623027699891626*softplus(0.9239332442221326 + 0.41157080299867177*$(x[1]) + 0.6856618842893334*$(x[2]) + 0.23250449051531152*$(x[3]) + -0.8171094215290706*$(x[4])) + 0.1823838704363938*softplus(-0.1441575648301785 + -0.7324773919132506*$(x[1]) + 0.1666248349273194*$(x[2]) + 0.5034013066229974*$(x[3]) + -0.3351866390106042*$(x[4]))))) - $q <= 0.0))

                     @objective(m, Min, q)

                     return m

max_stars_repo_name: PSORLab/RSActivationFunctions

max_issues_repo_name: PSORLab/RSActivationFunctions

max_forks_repo_name: PSORLab/RSActivationFunctions


text:
"""
Module Samsara:
-> The world for the agent: Day-to-day details and immediate focus of the agent
    -- The environment HUGIN is operating in Samsara: World mechanism!

Currently, this is a simple framework for how mechanics can be simulated:
    System mechanics is simulated by the function step_system_mechanics!(sys) som
    skriver til state (_observable_parameters or _latent_variables) til struct [sys]
    -> (1) Du kan overskrive step_system_mechanics! fra utsida:
    Dette er måten å lage ulike simuleringer.
    -> (2) Ved å deklarere fleire mutable struct SystemX <: AbstractSystem ∀ X,
    kan man også definere ulike system dynamikker ved overnevnte funksjon: step_system_mechanics(sys)
    -> (3) Struct har _observable_parameters og _latent_variables.
    Begge er tilgjengelig for ovenenvte funksjon!

    BAM! Skitbra!
"""
module Samsara

abstract type AbstractSystem end

#### Functions ###
" system_state(sys::AbstractSystem) returns system state of sys "
function system_state(sys::AbstractSystem)
    sys._observable_parameters  # return
end

" dimentionality(sys::AbstractSystem) returns dimentionality "
function dimensionality(sys::AbstractSystem)
    if isnothing(system_state(sys)[1]) # Ta bort?
        return 0
    else
        length(system_state(sys))
    end
end

include("demosys.jl")
#===========================================================================#
# Når man definerer eit nytt system, i ei anna fil, må denne koden takast
# inn i modulen i tekstform.
# -> No tar eg den med i module Samsara, siden denne skal brukes til TDD,
# men seinare kan du gjøre dette samme for nye system, lokalt i prosjektet.
# (include med local path, men siden den arver fra AbstractSystem, kan alt
# av Samsara's funksjoner brukes på den)
#===========================================================================#
include("circularsys.jl")
#===========================================================================#

export Samsara, AbstractSystem, dimensionality, system_state, step_system_mechanics!, set_action_in!           #EXPORT
#export EuclideanNode

include("linked_node.jl")
export LinkedCardinalNode

include("linked_gate.jl")
export LinkedGate

" convenience: Base.show for Pairs "
#function Base.show(io::IO, arg::Pair)
#    print(io, "Pair[ "*string(arg[1])*"  |||   "*string(arg[2])*" ]")
#end

end # module

max_stars_repo_name: leikanger/Samsara

max_issues_repo_name: leikanger/Samsara

max_forks_repo_name: leikanger/Samsara


text:
using AnalyticalEOR, Plots

swr = 0.2
sor = 0.2

kro0 = 1.0
krw0 = 0.5
nw = 3.0
no = 1.0

kr = RelPerms(swr, sor, krw0, kro0, nw, no)

μw = 1.0
μo = 3.0

si = 0.20
sj = 0.80

wf = solve_waterflooding(si, sj, kr, μw, μo)

tracer = solve_tracer(wf)

plot_tracer_fw(wf, tracer)

plot_sw_profile(wf, tracer, 0.1)
animate_sw_profile(wf, tracer)

max_stars_repo_path: examples/Ex2 - Tracers.jl

max_issues_repo_path: examples/Ex2 - Tracers.jl

max_forks_repo_path: examples/Ex2 - Tracers.jl


text:
#=
  Seven scientists in Turing.jl

  Port of the WebPPL model from
  https://mhtess.github.io/bdappl/chapters/03-simpleModels.html
  """
  Seven scientists with wildly-differing experimental skills all make a
  measurement of the same quantity. They get the answers

     x = {−27.020, 3.570, 8.191, 9.898, 9.603, 9.945, 10.056}.

  Intuitively, it seems clear that the first two scientists are pretty inept
  measurers, and that the true value of the quantity is probably just a bit
  below 10. The main problem is to find the posterior distribution over
  the measured quantity, telling us what we can infer from the measurement.
  A secondary problem is to infer something about the measurement skills of
  the seven scientists.

  """

  Using NUTS(1000,0.65):
  """
  Summary Statistics
  parameters      mean       std   naive_se      mcse         ess      rhat   ess_per_sec
      Symbol   Float64   Float64    Float64   Float64     Float64   Float64       Float64

          mu    9.1691    1.5779     0.0158    0.0388   1514.9796    1.0003      361.9158
   sigmas[1]   17.1826    2.1492     0.0215    0.0481   1812.3976    1.0024      432.9665
   sigmas[2]   10.1038    4.8808     0.0488    0.0947   2207.2242    1.0000      527.2872
   sigmas[3]    6.9711    5.3999     0.0540    0.1307   1840.5852    1.0000      439.7002
   sigmas[4]    5.7966    5.5111     0.0551    0.1452   1247.6705    1.0002      298.0579
   sigmas[5]    5.5648    5.4101     0.0541    0.1455   1266.5753    1.0008      302.5741
   sigmas[6]    5.6467    5.3356     0.0534    0.1272   1428.5790    1.0002      341.2754
   sigmas[7]    5.8013    5.4090     0.0541    0.1533   1491.1203    1.0003      356.2160

  Quantiles
  parameters      2.5%     25.0%     50.0%     75.0%     97.5%
      Symbol   Float64   Float64   Float64   Float64   Float64

          mu    4.4711    8.7121    9.6185    9.9814   11.3731
   sigmas[1]   12.1421   15.8454   17.5997   18.9317   19.8960
   sigmas[2]    2.3154    6.0300    9.5638   14.0701   19.2536
   sigmas[3]    0.4744    2.4507    5.4174   10.7118   18.7471
   sigmas[4]    0.1709    1.2017    3.7592    9.4112   18.4969
   sigmas[5]    0.1470    1.1639    3.5180    8.8308   18.5041
   sigmas[6]    0.1892    1.3070    3.7946    8.5708   18.6025
   sigmas[7]    0.1556    1.2915    3.7926    9.2773   18.1375
   """

=#

using Turing
include("jl_utils.jl")

@model function seven_scientists(data)
    n = length(data)
    mu ~ Normal(0, 30)
    sigmas ~ filldist(Uniform(0,20),n)
    for i in 1:n
        data[i] ~ Normal(mu, sigmas[i])
    end
end
data = [-27.020, 3.570, 8.191, 9.898, 9.603, 9.945, 10.056]

model = seven_scientists(data)

# chns = sample(model, Prior(), 10_000)
# chns = sample(model, MH(), 10_000)
# chns = sample(model, PG(15), 10_000)
# chns = sample(model, IS(), 10_000)
# chns = sample(model, SMC(), 10_000)
# chns = sample(model, SMC(), MCMCThreads(), 10_000, 4)

chns = sample(model, NUTS(1000,0.65), 10_000)
# chns = sample(model, HMC(0.1,10), 10_000)

display(chns)

max_stars_repo_name: tias/hakank

max_issues_repo_name: tias/hakank

max_forks_repo_name: tias/hakank


text:
module QuantumGate

    export σ⁰, σˣ, σʸ, σᶻ, P⁰, P¹, Had, Smat, Tmat, Rˣmat, Rʸmat, Rᶻmat
    include("pauli_matrices.jl")

    export ⊗
    include("kronecker_product.jl")

    export Xgate, Ygate, Zgate, Hadamardgate, Sgate, Tgate, Rˣgate, Rʸgate, Rᶻgate
    export controlled_Ugate, CNOTgate, SWAPgate
    include("gates.jl")

end


text:
function randcoords(nx, ny)
  K = zeros(Float32, 2, nx, ny)
  for j in 1:ny, i in 1:nx
    K[1,i,j] = rand() - 0.5
    K[2,i,j] = rand() - 0.5
  end
  K
end

function radialcoords(nr, narms; golden=false, offset=0)
  K = zeros(Float32, 2, nr, narms)
  for j in 1:narms
    θ = (j - 1 + offset) * π / (golden ? φ : narms)
    for i in 1:nr
      r = (i-1) / nr - 0.5
      K[1,i,j] = r*cos(θ)
      K[2,i,j] = r*sin(θ)
    end
  end
  K
end
function radialcoords(nr, pct::Float32=100.0; kwargs...)
  # can call with a radial percentage as well
  narms = iround(nr*pct / 100.0)
  radialcoords(nr, narms; kwargs...)
end

function radialcoords(nr, nprof, nleaves; pct=100.0, golden=false,
                      fid=false, roalt=false)
  nprof = iround(nprof*sqrt(pct/100.0))
  nleaves = iround(nleaves*sqrt(pct/100.0))
  K = zeros(Float32, 3, nr, nprof, nleaves)
  if golden
    α = 0.4656
    β = 0.6823
    ϕ = 0.0
    μ = 0.0
    for k in 1:nleaves, j in 1:nprof
      st = sqrt(1.0 - μ*μ)
      for i in 1:nr
        r = (i-1) / nr - 0.5
        K[1,i,j,k] = r*cos(ϕ)*st
        K[2,i,j,k] = r*sin(ϕ)*st
        K[3,i,j,k] = r*μ
      end
      ϕ = mod(ϕ + 2π*β, 2π)
      μ += 2α
      μ = μ > 1.0 ? μ - 2.0 : μ
    end
  else
    for j in 1:nprof
      if fid
        z = (2j - 1) / float(nprof)  - 1.0 # not sure about this
	z = j % 2 ? z : -z
      else
        z = (j - 0.5) / float(nprof) - 1.0
      end
      for k in 1:nleaves
        ϕ = sqrt(nprof * π / nleaves) * asin(z) + (k-1) * 2π / nleaves
	#from mmiffe_mxg.c phi = sqrt( `MP_radial_nr_angles * PI / `MP_radial_3d_koosh_interleaves ) * asin( z ) + interleaf * 2.0 * PI / `MP_radial_3d_koosh_interleaves;
        st = sqrt(1.0 - z*z)
        for i in 1:nr
          r = (i-1) / nr - 0.5
          K[1,i,j,k] = r*cos(ϕ)*st
          K[2,i,j,k] = r*sin(ϕ)*st
          K[3,i,j,k] = r*z
        end
      end
    end
  end
  reshape(K, 3, nr, nprof*nleaves)
end

function radialsdc(nr, nt; golden::Bool=false, offset::Int=0)
  # First guess radial sample density compensation
  # Uses Ram-Lak.  If golden angle, then the Ram-Lak
  # ramp is adjusted by the average neighbor distance.
  # sort angles in increasing order
  θ = (π / φ)*[offset:nt-1+offset] % (2π)
  Δθavg = 4π / nt
  θhat = sort(θ)
  isort = sortperm(θ)
  # find angular distance to neighbors
  Δθ = zeros(nt)
  Δθ[1] = θhat[2] - θhat[end] + 2π
  for t in 2:nt-1
    Δθ[t] = θhat[t+1] - θhat[t-1]
  end
  Δθ[end] = θhat[1] - θhat[end-1] + 2π
  weights = similar(Δθ)
  weights[isort] = Δθ / Δθavg
  r = linspace(1.0, 1.0/nr, int(nr/2))
  ramp = vcat(r, r[end:-1:1])
  if golden
    sdc = broadcast(*, ramp, weights')
  else
    sdc = repmat(ramp, 1, nt)
  end
  float32(sdc)
end

max_stars_repo_path: code/Julia Recon Scripts/coords.jl

max_issues_repo_path: code/Julia Recon Scripts/coords.jl

max_forks_repo_path: code/Julia Recon Scripts/coords.jl


text:
# Copyright (c) 2021, Eric Sabo
# All rights reserved.
#
# This source code is licensed under the BSD-style license found in the
# LICENSE file in the root directory of this source tree.

include("quantumcode.jl")
# include("tricolorcodes488trellis.jl")
# include("tricolorcodes666trellis.jl")

#############################
        # Misc codes
#############################

# miscellaneous small codes useful for testing functions

function fivequbitcode()
    # is a perfect code
    return QuantumCode(["XZZXI", "IXZZX", "XIXZZ", "ZXIXZ"])
end

function Q513()
    return fivequbitcode()
end

# should also do a test for other CSS construction via Hamming code and actually make that one default
function Steanecode()
    # standard
    return CSSCode(["XXXXIII", "XXIIXXI", "XIXIXIX", "ZZZZIII", "ZZIIZZI", "ZIZIZIZ"])

    # trellis numbering
    # return CSSCode(["XXIXXII", "IXXIXXI", "IIIXXXX", "ZZIZZII", "IZZIZZI", "IIIZZZZ"])
    # also ZZIZZII, ZIZZIZI, IZZZIIZ, XXIXXII, XIXXIXI, IXXXIIX
end

function Shorcode()
    return CSSCode(["ZZIIIIIII", "IZZIIIIII", "IIIZZIIII", "IIIIZZIII", "IIIIIIZZI", "IIIIIIIZZ", "XXXXXXIII", "IIIXXXXXX"])
end

function Q412()
    return CSSCode(["XXXX", "ZZII", "IIZZ"])
end

function Q422()
    return CSSCode(["XXXX", "ZZZZ"])
end

function Q511()
    return QuantumCode(["ZXIII", "XZXII", "IXZXI", "IIXZX"])
end

function Q823()
    F, _ = FiniteField(2, 1, "α")
    S = matrix(F, [1 0 0 0 1 0 0 0 1 1 1 1 0 0 0 0;
    0 0 0 1 0 1 0 0 1 0 0 0 0 1 0 0;
    0 1 0 0 1 1 1 0 0 0 1 1 1 0 1 0;
    0 0 1 0 1 1 1 0 0 1 1 0 1 1 0 0;
    0 0 1 1 1 0 1 0 0 0 0 1 0 1 1 1;
    0 0 0 0 0 0 1 1 0 0 1 0 0 0 1 0]);
    return QuantumCode(S, true)
end

function Q15RM()
    return QuantumCode(["ZIZIZIZIZIZIZIZ", "IZZIIZZIIZZIIZZ", "IIIZZZZIIIIZZZZ",
        "IIIIIIIZZZZZZZZ", "IIZIIIZIIIZIIIZ", "IIIIZIZIIIIIZIZ", "IIIIIZZIIIIIIZZ",
        "IIIIIIIIIZZIIZZ", "IIIIIIIIIIIZZZZ", "IIIIIIIIZIZIZIZ", "XIXIXIXIXIXIXIX",
        "IXXIIXXIIXXIIXX", "IIIXXXXIIIIXXXX", "IIIIIIIXXXXXXXX"])
end

#############################
#  Triangular Surface Codes #
#############################

function _triangularlattice(L::Int)
    # 0 - vertical
    # 1 - horizontal
    # 2 - diagonal
    numbering = zeros(Int, L, L, 3)
    num = 1
    for i in 1:L
        for j in 1:L
            for k in 1:3
                numbering[i, j, k] = num
                num += 1
            end
        end
    end
    return numbering
end

function _triangularlatticeXstabilizers(L::Int, numbering::Array{Int64, 3}, symp::Bool=true)
    F, _ = FiniteField(2, 1, "α")
    stabilizers = zero_matrix(F, L^2, 3 * L^2)
    r = 1
    for i in 1:L
        for j in 1:L
            stabilizers[r, numbering[i, j, 1]] = 1
            stabilizers[r, numbering[i, j, 2]] = 1
            stabilizers[r, numbering[i, j, 3]] = 1
            if i == 1
                stabilizers[r, numbering[end, j, 1]] = 1
            else
                stabilizers[r, numbering[i - 1, j, 1]] = 1
            end
            if j == 1
                stabilizers[r, numbering[i, end, 2]] = 1
            else
                stabilizers[r, numbering[i, j - 1, 2]] = 1
            end
            if i == 1 && j == 1
                stabilizers[r, numbering[end, end, 3]] = 1
            elseif i != 1 && j == 1
                stabilizers[r, numbering[i - 1, end, 3]] = 1
            elseif i == 1 && j != 1
                stabilizers[r, numbering[end, j - 1, 3]] = 1
            else
                stabilizers[r, numbering[i - 1, j - 1, 3]] = 1
            end
            r += 1
        end
    end
    if symp
        return hcat(stabilizers, zero_matrix(F, L^2, 3 * L^2))
    end
    return stabilizers
end

function _triangularlatticeZstabilizers(L::Int, numbering::Array{Int64, 3}, symp::Bool=true)
    F, _ = FiniteField(2, 1, "α")
    stabilizers = zero_matrix(F, 2 * L^2, 3 * L^2)
    r = 1
    for i in 1:L
        for j in 1:L
            stabilizers[r, numbering[i, j, 2]] = 1
            stabilizers[r, numbering[i, j, 3]] = 1
            if j == L
                stabilizers[r, numbering[i, 1, 1]] = 1
            else
                stabilizers[r, numbering[i, j + 1, 1]] = 1
            end

            r += 1
            stabilizers[r, numbering[i, j, 1]] = 1
            stabilizers[r, numbering[i, j, 3]] = 1
            if i == L
                stabilizers[r, numbering[1, j, 2]] = 1
            else
                stabilizers[r, numbering[i + 1, j, 2]] = 1
            end
            r += 1
        end
    end
    if symp
        return hcat(zero_matrix(F, 2 * L^2, 3 * L^2), stabilizers)
    end
    return stabilizers
end

function _triangularlatticeXlogicals(L::Int, numbering::Array{Int64, 3}, symp::Bool=true)
    # should be 0110110110
    z = zeros(UInt8, 3 * L^2)
    logical1 = zeros(UInt8, 3 * L^2)
    for j in 0:L - 1
        for k in 0:2
            if k == 0
                logical1[3 * j + k + 1] = 0x01
            elseif k == 1
                logical1[3 * j + k + 1] = 0x00
            else
                logical1[3 * j + k + 1] = 0x01
            end
        end
    end
    if symp
        logical1 = [logical1; z]
    end

    logical2 = zeros(UInt8, 3 * L^2)
    for j in 1:L
        for k in 0:2
            if k == 0
                logical2[3 * j + k + 1] = 0x01
            elseif k == 1
                logical2[3 * j + k + 1] = 0x00
            else
                logical2[3 * j + k + 1] = 0x01
            end
        end
    end
    if symp
        logical2 = [logical2; z]
    end

  return [logical1, logical2]
end

function _triangularlatticeZlogicals(L::Int, numbering::Array{Int64, 3}, symp::Bool=true)
    # should be 1001001001
    x = zeros(UInt8, 3 * L^2)
    logical1 = zeros(UInt8, 3 * L^2)
    for j in 0:L - 1
        for k in 0:2
            if k == 0
                logical1[3 * j + k + 1] = 0x01
            elseif k == 1
                logical1[3 * j + k + 1] = 0x00
            else
                logical1[3 * j + k + 1] = 0x01
            end
        end
    end
    if symp
        logical1 = [z; logical1]
    end

    logical2 = zeros(UInt8, 3 * L^2)
    for j in 1:L
        for k in 0:2
            if k == 0
                logical2[3 * j + k + 1] = 0x01
            elseif k == 1
                logical2[3 * j + k + 1] = 0x00
            else
                logical2[3 * j + k + 1] = 0x01
            end
        end
    end
    if symp
        logical2 = [z; logical2]
    end

  return [logical1, logical2]
end

function triangularsurfacecode(L::Int)
    numbering = _triangularlattice(L)
    Xstabs = _triangularlatticeXstabilizers(L, numbering, false)
    # println(rank(Xstabs))
    Zstabs = _triangularlatticeZstabilizers(L, numbering, false)
    # println(Zstabs)
    # logicals = [triangularlatticeXlogicals(L, numbering), triangularlatticeZlogicals(L, numbering)]
    return CSSCode(Xstabs[1:end - 1, :], Zstabs[1:end - 1, :])
end

#############################
#  Rotated Surface Codes #
#############################

function _RSurfstabslogs(d::Int)
    n = d^2
    E, ω = FiniteField(2, 2, "ω")
    X = E(1)
    S = zero_matrix(E, n - 1, n)
    row = 1

    # X's
    i = 1
    while i <= n - d
        S[row, i] = X
        S[row, i + 1] = X
        S[row, i + d] = X
        S[row, i + d + 1] = X
        row += 1
        if (i + 2) % d == 0
            i += 4
        else
            i += 2
        end
    end

    # top row X's
    i = 2
    while i <= d - 1
        S[row, i] = X
        S[row, i + 1] = X
        row += 1
        i += 2
    end

    # bottom row X's
    i = d * (d - 1) + 1
    while i <= d * d - 2
        S[row, i] = X
        S[row, i + 1] = X
        row += 1
        i += 2
    end

    # Z's
    i = 2
    while i < n - d
        S[row, i] = ω
        S[row, i + 1] = ω
        S[row, i + d] = ω
        S[row, i + d + 1] = ω
        row += 1
        if (i + 2) % d == 0
            i += 4
        else
            i += 2
        end
    end

    # left Z's
    i = 1
    while i < d * (d - 1)
        S[row, i] = ω
        S[row, i + d] = ω
        row += 1
        i += 2 * d
    end

    # right Z's
    i = 2 * d
    while i < d * d
        S[row, i] = ω
        S[row, i + d] = ω
        row += 1
        i += 2 * d
    end

    logs = zero_matrix(E, 2, n)
    i = d
    while i <= d * d
        logs[1, i] = 0x01
        i += d
    end
    i = 1
    while i <= d
        logs[2, i] = ω
        i += 1
    end

    return S, logs
end

# flint is being ridiculous here
# julia> @time Q = rotatedsurfacecode(21);
# 226.128107 seconds (1.73 G allocations: 80.659 GiB, 18.33% gc time, 0.09% compilation time)
# julia> Base.summarysize(Q)
# 16872
function rotatedsurfacecode(d::Int)
    d >= 3 || error("Current implementation requires d ≥ 3.")

    S, logs = _RSurfstabslogs(d)
    Q = QuantumCode(S)
    setlogicals!(Q, logs)
    return Q
end

#############################
#     XZZX Surface Codes    #
#############################

function _XZZXstabslogs(d::Int)
    n = d^2
    E, ω = FiniteField(2, 2, "ω")
    S = zero_matrix(E, n - 1, n)
    row = 1
    X = E(1)

    i = 1
    for i in 1:n - d
        if i % d != 0
            S[row, i] = X
            S[row, i + 1] = ω
            S[row, i + d] = ω
            S[row, i + d + 1] = X
            row += 1;
        end
    end

    # top row ZX's
    i = 2
    while i <= d - 1
        S[row, i] = ω
        S[row, i + 1] = X
        row += 1
        i += 2
    end

    # bottom row XZ's
    i = d * (d - 1) + 1
    while i <= d * d - 2
        S[row, i] = X
        S[row, i + 1] = ω
        row += 1
        i += 2
    end

    # left ZX's
    i = 1
    while i < d * (d - 1)
        S[row, i] = ω
        S[row, i + d] = X
        row += 1
        i += 2 * d
    end

    # right XZ's
    i = 2 * d
    while i < d * d
        S[row, i] = X
        S[row, i + d] = ω
        row += 1
        i += 2 * d
    end

    logs = zero_matrix(E, 2, n)
    i = d
    count = 1
    while i <= d * d
        if count % 2 == 1
            logs[1, i] = X
        else
            log[1, i] = ω
        end
        i += d
        count += 1
    end
    i = 1
    count = 1
    while i <= d
        if count % 2 == 1
            log[2, i] = ω
        else
            log[2, i] = X
        end
        i += 1
        count += 1
    end

    return S, logs
end

function XZZXsurfacecode(d::Int)
    d >= 3 || error("Current implementation requires d ≥ 3.")

    S, logs = _XZZXstabslogs(d)
    Q = QuantumCode(S)
    setlogicals!(Q, logs)
    return Q
end

################################
# Triangular Color Codes 4.8.8 #
################################

function tricolorcode488(d::Int)
    3 <= d <= 21 || error("Current implementation requires 3 ≤ d ≤ 21.")

    if d == 3
        S, logs = _488d3trellis()
        Q = QuantumCode(S)
        setlogicals!(Q, logs)
        return Q
    elseif d == 5
        S, logs = _488d5trellis()
        Q = QuantumCode(S)
        setlogicals!(Q, logs)
        return Q
    elseif d == 7
        S, logs = _488d7trellis()
        Q = QuantumCode(S)
        setlogicals!(Q, logs)
        return Q
    elseif d == 9
        S, logs = _488d9trellis()
        Q = QuantumCode(S)
        setlogicals!(Q, logs)
        return Q
    elseif d == 11
        S, logs = _488d11trellis()
        Q = QuantumCode(S)
        setlogicals!(Q, logs)
        return Q
    elseif d == 13
        S, logs = _488d13trellis()
        Q = QuantumCode(S)
        setlogicals!(Q, logs)
        return Q
    elseif d == 15
        S, logs = _488d15trellis()
        Q = QuantumCode(S)
        setlogicals!(Q, logs)
        return Q
    elseif d == 17
        S, logs = _488d17trellis()
        Q = QuantumCode(S)
        setlogicals!(Q, logs)
        return Q
    elseif d == 19
        S, logs = _488d19trellis()
        Q = QuantumCode(S)
        setlogicals!(Q, logs)
        return Q
    elseif d == 21
        S, logs = _488d21trellis()
        Q = QuantumCode(S)
        setlogicals!(Q, logs)
        return Q
    end
end

################################
# Triangular Color Codes 6.6.6 #
################################

function tricolorcode666(d::Int)
    3 <= d <= 21 || error("Current implementation requires 3 ≤ d ≤ 21.")

    if d == 3
        # same as 4.8.8
        S, logs = _488d3trellis()
        Q = QuantumCode(S)
        setlogicals!(Q, logs)
        return Q
    elseif d == 5
        S, logs = _666d5trellis()
        Q = QuantumCode(S)
        setlogicals!(Q, logs)
        return Q
    elseif d == 7
        S, logs = _666d7trellis()
        Q = QuantumCode(S)
        setlogicals!(Q, logs)
        return Q
    elseif d == 9
        S, logs = _666d9trellis()
        Q = QuantumCode(S)
        setlogicals!(Q, logs)
        return Q
    elseif d == 11
        S, logs = _666d11trellis()
        Q = QuantumCode(S)
        setlogicals!(Q, logs)
        return Q
    elseif d == 13
        S, logs = _666d13trellis()
        Q = QuantumCode(S)
        setlogicals!(Q, logs)
        return Q
    elseif d == 15
        S, logs = _666d15trellis()
        Q = QuantumCode(S)
        setlogicals!(Q, logs)
        return Q
    elseif d == 17
        S, logs = _666d17trellis()
        Q = QuantumCode(S)
        setlogicals!(Q, logs)
        return Q
    elseif d == 19
        S, logs = _666d19trellis()
        Q = QuantumCode(S)
        setlogicals!(Q, logs)
        return Q
    elseif d == 21
        S, logs = _666d21trellis()
        Q = QuantumCode(S)
        setlogicals!(Q, logs)
        return Q
    end
end

max_stars_repo_name: esabo/CodingTheory

max_issues_repo_name: esabo/CodingTheory

max_forks_repo_name: esabo/CodingTheory


text:
@testset begin "Connectivities"
  # basic tests
  c = connect((1,2,3), Triangle)
  @test issimplex(c)
  @test paramdim(c) == 2
  @test indices(c) == (1,2,3)
  @test materialize(c, P2[(0,0),(1,0),(0,1)]) == Triangle(P2(0,0), P2(1,0), P2(0,1))

  # tuple from other collections
  c = connect(Tuple([1,2,3]), Triangle)
  @test paramdim(c) == 2
  @test indices(c) == (1,2,3)
  @test materialize(c, P2[(0,0),(1,0),(0,1)]) == Triangle(P2(0,0), P2(1,0), P2(0,1))

  # incorrect number of vertices for polytope
  @test_throws AssertionError connect((1,2,3,4), Triangle)
  @test_throws AssertionError connect((1,2,3), Quadrangle)

  # heterogeneous collections
  c = connect.([(1,2,6,5),(2,4,6),(4,3,5,6),(1,5,3)], Ngon)
  @test c[1] isa Connectivity{Quadrangle}
  @test c[2] isa Connectivity{Triangle}
  @test c[3] isa Connectivity{Quadrangle}
  @test c[4] isa Connectivity{Triangle}
  @test issimplex.(c) == [false,true,false,true]

  # ommitting polytope type means polygon or segment
  @test connect((1,2)) isa Connectivity{Segment}
  @test connect((1,2,3)) isa Connectivity{Triangle}
  @test connect((1,2,3,4)) isa Connectivity{Quadrangle}
  @test connect((1,2,3,4,5)) isa Connectivity{Pentagon}
  @test connect((1,2,3,4,5,6)) isa Connectivity{Hexagon}

  # polyhedron connectivities
  c = connect((1,2,3,4), Tetrahedron)
  @test issimplex(c)
  @test paramdim(c) == 3
  @test indices(c) == (1,2,3,4)
  points = P3[(0,0,0),(1,0,0),(0,1,0),(0,0,1)]
  @test materialize(c, points) == Tetrahedron(points...)
end


text:
using Requires
@init @require CUDAnative = "be33ccc6-a3ff-5ff2-a52e-74243cff1e17" begin
  using .CUDAnative
end

function update!(Q, offset, ::Val{i}, yn, ΔYnj, fYnj, αi, βi, γi, d_i, dt) where i
  @inbounds @loop for e = (1:length(Q);
                           (blockIdx().x - 1) * blockDim().x + threadIdx().x)
    if i > 2
      ΔYnj[i-2][e] = Q[e] - yn[e] # is 0 for i == 2
    end
    Q[e] = yn[e] # (1a)
    offset[e] = (βi[1]/d_i) .* fYnj[1][e] # (1b)
    @unroll for j = 2:i-1
      Q[e] += αi[j] .* ΔYnj[j-1][e] # (1a cont.)
      offset[e] += (γi[j]/(d_i*dt)) * ΔYnj[j-1][e] + (βi[j]/d_i) * fYnj[j][e] # (1b cont.)
    end
  end
end

max_stars_repo_name: akshaysridhar/CLIMA

max_issues_repo_name: akshaysridhar/CLIMA

max_forks_repo_name: akshaysridhar/CLIMA


text:
#using PyCall;
#using PyPlot;

module mnist

export read_mnist_data, show_some_imgs

# Read images from file
function read_mnist_img_file(path::String, correct_magic_number::Int32)
    img_file = open(path, "r")
    magic_number = bswap(read(img_file, Int32));
    @assert magic_number == correct_magic_number
    n = int64(bswap(read(img_file, Int32)));
    rows = int64(bswap(read(img_file, Int32)));
    columns = int64(bswap(read(img_file, Int32)));
    data = float64(read(img_file, Uint8, (rows * columns, n))) ./ 255.0
    close(img_file);
    return data;
end

# Read labels
function read_mnist_label_file(path::String, correct_magic_number::Int32)
    label_file = open(path, "r")
    magic_number = bswap(read(label_file, Int32));
    @assert magic_number == correct_magic_number
    n = int64(bswap(read(label_file, Int32)));
    data = read(label_file, Uint8, (n))
    targets = zeros(10, n)
    for i = 1:n
        targets[data[i]+1, i] = 1
    end
    close(label_file);
    return targets;
end

function read_mnist_data()
    # Read images
    train_img = read_mnist_img_file("data/train-images.idx3-ubyte", int32(2051))
    test_img = read_mnist_img_file("data/t10k-images.idx3-ubyte", int32(2051))
    # Read labels
    train_label = read_mnist_label_file("data/train-labels.idx1-ubyte",
                                        int32(2049))
    test_label = read_mnist_label_file("data/t10k-labels.idx1-ubyte",
                                       int32(2049))
    return ["train_img" => train_img,
            "test_img" => test_img,
            "train_label" => train_label,
            "test_label" => test_label]
end

function show_some_imgs(imgs, labels, rows::Integer, cols::Integer)
    d = zeros(rows, cols)

    index1 = ceil(rand()*size(imgs,2))
    row = permutedims(reshape(imgs[:,index1], 28, 28), [2,1])
    d[1,1] = find(labels[:,index1])[1]
    for j = 2:cols
        index1 = ceil(rand()*size(imgs,2))
        row = hcat(row, permutedims(reshape(imgs[:,index1], 28, 28), [2,1]))
        d[1,j] = find(labels[:,index1])[1]
    end

    m = row

    for i = 2:rows
        index1 = ceil(rand()*size(imgs,2))
        row = permutedims(reshape(imgs[:,index1], 28, 28), [2,1])
        d[i,1] = find(labels[:,index1])[1]
        for j = 2:cols
            index1 = ceil(rand()*size(imgs,2))
            row = hcat(row, permutedims(reshape(imgs[:,index1], 28, 28), [2,1]))
            d[i,j] = find(labels[:,index1])[1]
        end
        m = vcat(m, row)
    end
    println(int(d .- 1))
    imshow(m)
end

end

max_stars_repo_name: tudor-berariu/ml-hacks

max_issues_repo_name: tudor-berariu/ml-hacks

max_forks_repo_name: tudor-berariu/ml-hacks