text:
using Distributions

# simulates the shell compression process, based on the number of sorted out samples in the ns iterations
# result is relative to volume of original polytope
function sim_shell_compression( itd::Vector{NSIterData} ; repetitions::Int64=-1)
      # default value
      if(repetitions<0); repetitions = 1000; end

      n_a = map( xi -> xi.n_live_samples_a , itd)
      n_b = map( xi -> xi.n_live_samples_b , itd)

      log_xv      = NaN * ones(repetitions,length(n_a)+1)
      log_xv[:,1] = 0.0

      x_betas = NaN * ones(repetitions,length(n_a))

      # draw betas..
      for zi in 1:length(n_a)

            # k = sorted out, n = samples before sorting out
            # then shrinking is  1-beta( k , n+1-k )
            x_betas[:,zi]  = ( 1.0 - rand(Distributions.Beta( n_a[zi]-n_b[zi] , n_b[zi]+1) , repetitions ) )
            log_xv[:,zi+1] = log_xv[:,zi] + log.(x_betas[:,zi])
      end
      x_v    = exp.(log_xv)
      x_mean = mean( x_v , 1 )
      return ( x_mean , x_v , log_xv )
end

"""
  eval_integral( volume::Float64, itd::Vector{NSIterData} ; repetitions::Int64=-1)

returns (Float64,Vector{Float64},Array{Float64,2})
where [1] is the estimated mean, [2] is the sampled integration results
and [3] is the (sampled) integrated parts of the ns iterations

evaluates the integral based on NS data, for the pdf with pdf(x) = exp( - fcost(x))
volume is the volume of the sampled space
"""
function eval_integral( volume::Float64, itd::Vector{NSIterData} ; repetitions::Int64=-1)
      ssc_1 = sim_shell_compression(itd;repetitions=repetitions)

      ssc = volume * ssc_1[2]

      # consider sorted out ns samples as mc sample of the shell
      iv_parts = Array{Float64}(size(ssc,1),length(itd))
      for zi=1:length(itd)
            mc_xf        = mean( exp.( - itd[zi].Xf_sorted_out )  )
            iv_parts[:,zi] = ( ssc[:,zi]-ssc[:,zi+1] ) * mc_xf
      end

      iv_sum = sum(iv_parts,2)
      return (mean(iv_sum),iv_sum,iv_parts)
end

"""
  resample_ns_data( itd::Vector{NSIterData} )

returns Array{Float64,2}
draws a sample of the target distribution, given the output of a ns run
"""
function resample_ns_data( itd::Vector{NSIterData} )
      int_result = eval_integral( 1.0 , itd )
      iv_parts   = int_result[3]

      # now draw proportional to iv_parts from the sorted out samples of the iterations..
      X   = zeros(size(itd[1].X_sorted_out,1),0)
      X_p = zeros(0)
      for zi in 1:length(itd)
            X   = [ X     itd[zi].X_sorted_out ]
            X_p = [ X_p ; ones(size(itd[zi].X_sorted_out,2)) * (mean(iv_parts[:,zi]) / size(itd[zi].X_sorted_out,2) )  ]
      end

      # sample until we hit collision:
      sampled_indeces = Vector{Int64}()
      no_collision = true
      while(no_collision)
            si = SamplingHS3.sample_proportional( X_p , 4 )
            append!( sampled_indeces , si )
            no_collision = length(unique(sampled_indeces))==length(sampled_indeces)
      end
      return ( X[ : ,  sampled_indeces] , X , X_p  )
end


text:
using MakieLayout
using Makie

using AbstractPlotting: px

begin
    scene = Scene(resolution = (1000, 1000));
    screen = display(scene)
    campixel!(scene);

    nrows = 4
    ncols = 4

    maingl = GridLayout(scene, 1, 1, alignmode=Outside(30))

    gridgl = maingl[1, 1] = GridLayout(
        nrows, ncols)

    las = [gridgl[i, j] = LAxis(scene) for i in 1:nrows, j in 1:ncols]

    linkxaxes!(las...)
    linkyaxes!(las...)

    for i in 1:nrows, j in 1:ncols

        scatter!(las[i, j], rand(200, 2) .+ [i j], markersize=20px, color=(:black, 0.3))

        i > 1 && (las[i, j].titlevisible = false)
        j > 1 && (las[i, j].ylabelvisible = false)
        j > 1 && (las[i, j].yticklabelsvisible = false)
        j > 1 && (las[i, j].yticksvisible = false)
        i < nrows && (las[i, j].xticklabelsvisible = false)
        i < nrows && (las[i, j].xticksvisible = false)
        las[i, j].xlabelvisible = false
    end

    gridgl[:, end, Right()] = LText(scene, text="Side Protrusion", rotation=-pi/2, padding=(20, 0, 0, 0), halign=:right)
    gridgl[end, :, Bottom()] = LText(scene, text="x label", padding=(0, 0, 0, 50), valign=:bottom)

    maingl[2, 1] = LAxis(scene)
    maingl[:, 2] = LAxis(scene)

    gridgl[:, :, TopLeft()] = LText(scene, text="A", valign=:top, halign=:left, textsize=50, padding=(0, 30, 30, 0))
    maingl[2, 1, TopLeft()] = LText(scene, text="B", valign=:top, halign=:left, textsize=50, padding=(0, 30, 30, 0))
    maingl[:, 2, TopLeft()] = LText(scene, text="C", valign=:top, halign=:left, textsize=50, padding=(0, 30, 30, 0))

    tl = maingl[0, :] = LText(scene, text="Super Title", textsize=50)
    stl = maingl[2:end, end+1] = LText(scene, text="Side Title", textsize=50, rotation=-pi/2)

    slgl = maingl[end+1, 1:end-1] = GridLayout(1, 2)

    slgl[1, 1] = LText(scene, text="Supertitle Size", halign=:left)
    sl1 = slgl[1, 2] = LSlider(scene, height = 30, range = 1:200)
    on(sl1.value) do val
        tl.attributes.textsize = val
    end

    slgl[2, 1] = LText(scene, text="Sidetitle Size", halign=:left)
    sl2 = slgl[2, 2] = LSlider(scene, height = 30, range = 1:200)
    on(sl2.value) do val
        stl.attributes.textsize = val
    end

    nothing
end

ext: jl


text:
using Test

enabled_tests = lowercase.(ARGS)
function addtests(fname)
  key = lowercase(splitext(fname)[1])
  if isempty(enabled_tests) || key in enabled_tests
    include(fname)
  end
end

addtests("models/RandomForestClassifier.jl")
addtests("models/Stats.jl")
addtests("models/Histogram.jl")
addtests("models/GaussianMixture.jl")
addtests("models/KMeans.jl")
addtests("models/GMeans.jl")
addtests("models/HierachicalClustering.jl")
addtests("controllers/BasicStatisticsController.jl")
addtests("controllers/CorrelationsController.jl")
addtests("controllers/GaussianMixturesController.jl")
addtests("controllers/TimeSeriesController.jl")
addtests("controllers/KMeansController.jl")
addtests("controllers/GMeansController.jl")
addtests("controllers/HierachicalClusteringsController.jl")


text:
using Traceur, Test

warns_for(ws, x) = any(w -> occursin(x, w.message), ws)
warns_for(ws, x, xs...) = warns_for(ws, x) && warns_for(ws, xs...)

y = 1
const cy = 2
naive_relu(x) = x < 0 ? 0 : x

randone() = rand() < 0.5 ? Int(1) : rand() < 0.5 ? Float64(1) : BigInt(1)

function naive_sum(xs)
  s = 0
  for x in xs
    s += x
  end
  return s
end

f(x) = x+y

function f2(x)
  foo = y
  sin(x)+y
end

g(x) = x+cy

naive_sum_wrapper(x) = naive_sum(x)

x = 1
my_add(y) = x + y

module Foo
module Bar
function naive_sum(xs)
  s = 0
  for x in xs
    s += x
  end
  return s
end
end
naive_sum_wrapper(x) = Bar.naive_sum(x)
end

@should_not_warn my_stable_add(y) = my_add(y)
my_stable_add_undecorated(y) = my_add(y)

@testset "Traceur" begin
  ws = Traceur.warnings(() -> naive_relu(1))
  @test isempty(ws)

  ws = Traceur.warnings(() -> naive_relu(1.0))
  @test warns_for(ws, "returns")

  ws = Traceur.warnings(() -> randone())
  @test warns_for(ws, "returns")

  ws = Traceur.warnings(() -> naive_sum([1]))
  @test isempty(ws)

  ws = Traceur.warnings(() -> naive_sum([1.0]))
  @test warns_for(ws, "assigned", "returns")

  ws = Traceur.warnings(() -> f(1))
  @test warns_for(ws, "global", "dispatch", "returns")

  ws = Traceur.warnings(() -> f2(1))
  @test warns_for(ws, "global", "dispatch", "returns")

  ws = Traceur.warnings(() -> g(1))
  @test isempty(ws)

  @testset "depth limiting" begin
    ws = Traceur.warnings(() -> naive_sum_wrapper(rand(3)); maxdepth = 0)
    @test length(ws) == 1
    @test warns_for(ws, "returns")

    ws = Traceur.warnings(() -> naive_sum_wrapper(rand(3)); maxdepth = 1)
    @test length(ws) == 4
    @test warns_for(ws, "assigned", "returns")
  end

  @testset "module specific" begin
    ws = Traceur.warnings(() -> Foo.naive_sum_wrapper(rand(3)); maxdepth = 2, modules=[Foo])
    @test length(ws) == 1
    @test warns_for(ws, "returns")

    ws = Traceur.warnings(() -> Foo.naive_sum_wrapper(rand(3)); maxdepth = 2, modules=[Foo.Bar])
    @test length(ws) == 3
    @test warns_for(ws, "assigned", "returns")

    ws = Traceur.warnings(() -> Foo.naive_sum_wrapper(rand(3)); maxdepth = 2, modules=[Foo, Foo.Bar])
    @test length(ws) == 4
    @test warns_for(ws, "assigned", "returns")
  end

  @testset "test utilities" begin
    @test_nowarn @check my_add(1)
    @test_throws AssertionError @check my_stable_add(1)
    @test_throws AssertionError @check my_stable_add_undecorated(1) nowarn=[my_stable_add_undecorated]
  end
end


text:
using Distributions
using LaTeXStrings

# Create plot of Beta(α,β) density
α = 10
β = 5

x_values = 0:0.01:1
y_values = [pdf(Beta(α,β),x) for x in x_values]

plot(x_values, y_values, seriestype = :line, legend = false,
 yaxis = L"f(\theta)", xaxis = L"\theta", grid = false)

savefig("beta_prior_coin_ex.svg")

max_stars_repo_name: seanrattana/courses

max_issues_repo_name: seanrattana/courses

max_forks_repo_name: seanrattana/courses


text:
function pc_range(die)
    low = extract_attribute(die, DWARF.DW_AT_low_pc)
    high = extract_attribute(die, DWARF.DW_AT_high_pc)
    (isnull(low) || isnull(high)) && return Nullable{UnitRange{UInt64}}()
    low = get(low).value
    high = get(high)
    if high.spec.form == DW_FORM_addr
        high = high.value
    else
        high = low + high.value
    end
    Nullable{UnitRange{UInt64}}(low:high)
end

function _searchforip(it, ip)
    for x in it
        pcr = pc_range(x)
        isnull(pcr) && continue
        pcr = get(pcr)
        (first(pcr) <= ip < last(pcr)) && return x
    end
    error("Not found")
end

searchcuforip(dbgs, ip) =
  _searchforip(UnitIterator{DWARFCUHeader}(dbgs.debug_info, dbgs.debug_abbrev), ip)
searchspforip(cu, ip) = _searchforip(ChildIterator(cu), ip)


text:
import Base.Git
include("gitutils.jl")

dir = string("tmp.",randstring())
@test !ispath(dir)
mkdir(dir)
@test isdir(dir)
try cd(dir) do

    run(`git init -q`)
    run(`git commit -q --allow-empty -m "initial empty commit"`)
    git_verify(Dict(), Dict(), Dict())

    # each path can have one of these content in each of head, index, work
    # for a total of length(contents)^3 = 4^3 = 64 combinations.
    # each path can be in any of these 64 "superpositions" before & after
    # for a total of 64^2 = 4096 files needed to test all transitions
    # between before and after superpositions of git repo states.

    contents = [nothing, "foo", "bar", {"baz"=>"qux"}]
    b = length(contents)
    states = [ [ base(b,k,6) => contents[rem(div(k,b^p),b)+1] for k=0:(b^3)^2-1 ] for p=0:5 ]

    git_setup(states[1:3]...)
    try Git.transact() do
        git_setup(states[4:6]...)
        throw(nothing)
    end catch x
        is(x,nothing) || rethrow()
    end
    git_verify(states[1:3]...)

end finally run(`rm -rf $dir`) end

max_stars_repo_name: porterjamesj/julia

max_issues_repo_name: porterjamesj/julia

max_forks_repo_name: porterjamesj/julia


text:
function construct_network!(ps_m::CanonicalModel,
                            system_formulation::Type{CopperPlatePowerModel},
                            sys::PSY.System,
                            time_range::UnitRange{Int64}; kwargs...)

    bus_count = length(sys.buses)

    copper_plate(ps_m, :var_active, bus_count, time_range)

    return
end

function construct_network!(ps_m::CanonicalModel,
                            system_formulation::Type{StandardPTDFForm},
                            sys::PSY.System,
                            time_range::UnitRange{Int64}; kwargs...)

    if :PTDF in keys(kwargs)

        ac_branches = [br for br in sys.branches if !isa(br, PSY.DCLine)]

        # TODO: Get DC Lines and model DC Lines

        flow_variables(ps_m, system_formulation, ac_branches, time_range)

        ptdf_networkflow(ps_m, ac_branches, sys.buses, :var_active, kwargs[:PTDF], time_range)

    else
        throw(ArgumentError("no PTDF matrix supplied"))
    end

    return

end

function construct_network!(ps_m::CanonicalModel,
                            system_formulation::Type{S},
                            sys::PSY.System,
                            time_range::UnitRange{Int64}; kwargs...) where {S <: PM.AbstractPowerFormulation}

    powermodels_network!(ps_m, system_formulation, sys, time_range)

    return

end


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

module Git
#
# some utility functions for working with git repos
#
using Base: shell_escape
using VersionParsing
export gitcmd # determined by deps/build.jl and saved in deps/deps.jl

depsjl = joinpath(dirname(@__FILE__), "..", "deps", "deps.jl")
isfile(depsjl) ? include(depsjl) : error("Git.jl not properly installed. " *
    "Please run\nPkg.build(\"Git\")")

"""
    Git.dir(d)

Return the path to the default `.git` for the given repository directory, or the
path to use in place of the default `.git`.
"""
function dir(d)
    g = joinpath(d,".git")
    isdir(g) && return g
    normpath(d, Base.readchomp(setenv(`$gitcmd rev-parse --git-dir`, dir=d)))
end

"""
    Git.git([d])

Return a Git command that refers to the work tree and directory given by `d`, or the
current work tree and directory if `d` is not specified.
"""
git() = gitcmd
function git(d)
    isempty(d) && return gitcmd
    work_tree = abspath(d)
    git_dir = joinpath(work_tree, dir(work_tree))
    normpath(work_tree, ".") == normpath(git_dir, ".") ? # is it a bare repo?
    `$gitcmd --git-dir=$work_tree` : `$gitcmd --work-tree=$work_tree --git-dir=$git_dir`
end

"""
    Git.cmd(args; dir="")

Return a Git command from the given arguments, acting on the repository given in `dir`.
"""
cmd(args::Cmd; dir="") = `$(git(dir)) $args`

"""
    Git.run(args; dir="", out=stdout)

Execute the Git command from the given arguments `args` on the repository `dir`, writing
the results to the output stream `out`.
"""
run(args::Cmd; dir="", out=stdout) = Base.run(pipeline(cmd(args,dir=dir), out))

"""
    Git.readstring(args; dir="")

Read the result of the Git command using the given arguments on the given repository
as a string.
"""
readstring(args::Cmd; dir="") = read(cmd(args,dir=dir), String)

"""
    Git.readchomp(args; dir="")

Read the result of the Git command using the given arguments on the given repository
as a string, removing a single trailing newline if present.
"""
readchomp(args::Cmd; dir="") = Base.readchomp(cmd(args,dir=dir))

"""
    Git.success(args; dir="")

Determine whether the Git command using the given arguments on the given repository
executed successfully.
"""
function success(args::Cmd; dir="")
    g = git(dir)
    Base.readchomp(`$g rev-parse --is-bare-repository`) == "false" &&
        Base.run(`$g update-index -q --really-refresh`)
    Base.success(`$g $args`)
end

"""
    Git.version()

Return the version of Git being used by the package.
"""
version() = vparse(readchomp(`version`))

"""
    Git.modules(args; dir="")

Apply the Git command with the given arguments on the given repository to the
configuration file `.gitmodules` and read the result as a string.
"""
modules(args::Cmd; dir="") = readchomp(`config -f .gitmodules $args`, dir=dir)

"""
    Git.different(verA, verB, path; dir="")

Determine whether two trees are different with respect to the given path.
"""
different(verA::AbstractString, verB::AbstractString, path::AbstractString; dir="") =
    !success(`diff-tree --quiet $verA $verB -- $path`, dir=dir)

"""
    Git.dirty([paths]; dir="")

Determine whether the paths in the given repository are dirty, i.e. contain modified but
uncommitted tracked files.
"""
dirty(; dir="") = !success(`diff-index --quiet HEAD`, dir=dir)
dirty(paths; dir="") = !success(`diff-index --quiet HEAD -- $paths`, dir=dir)

"""
    Git.staged([paths]; dir="")

Determine whether the paths in the given repository contain staged files.
"""
staged(; dir="") = !success(`diff-index --quiet --cached HEAD`, dir=dir)
staged(paths; dir="") = !success(`diff-index --quiet --cached HEAD -- $paths`, dir=dir)

"""
    Git.unstaged([paths]; dir="")

Determine whether the paths in the given repository contain unstaged files.
"""
unstaged(; dir="") = !success(`diff-files --quiet`, dir=dir)
unstaged(paths; dir="") = !success(`diff-files --quiet -- $paths`, dir=dir)

"""
    Git.iscommit(name; dir="")

Determine whether `name` refers to a commit in the repository `dir`. `name` can be a
single SHA1 or a vector of SHA1s.
"""
iscommit(name; dir="") = success(`cat-file commit $name`, dir=dir)
function iscommit(sha1s::Vector; dir="")
    indexin(sha1s,split(readchomp(`log --all --format=%H`, dir=dir),"\n")).!=nothing
end

"""
    Git.attached(; dir="")

Determine whether HEAD is attached to a commit in the given respository.
"""
attached(; dir="") = success(`symbolic-ref -q HEAD`, dir=dir)

"""
    Git.branch(; dir="")

Return the name of the current active branch in the given repository.
"""
branch(; dir="") = readchomp(`rev-parse --symbolic-full-name --abbrev-ref HEAD`, dir=dir)

"""
    Git.head(; dir="")

Return the commit to which HEAD currently refers.
"""
head(; dir="") = readchomp(`rev-parse HEAD`, dir=dir)

struct State
    head::String
    index::String
    work::String
end

"""
    Git.snapshot(; dir="")

Return a `State` object that captures a snapshot of the given repository.
"""
function snapshot(; dir="")
    head = readchomp(`rev-parse HEAD`, dir=dir)
    index = readchomp(`write-tree`, dir=dir)
    work = try
        if length(readdir(abspath(dir))) > 1
            run(`add --all`, dir=dir)
            run(`add .`, dir=dir)
        end
        readchomp(`write-tree`, dir=dir)
    finally
        run(`read-tree $index`, dir=dir) # restore index
    end
    State(head, index, work)
end

"""
    Git.restore(s::State; dir="")

Restore the given repository to the state `s`.
"""
function restore(s::State; dir="")
    run(`reset -q --`, dir=dir)               # unstage everything
    run(`read-tree $(s.work)`, dir=dir)       # move work tree to index
    run(`checkout-index -fa`, dir=dir)        # check the index out to work
    run(`clean -qdf`, dir=dir)                # remove everything else
    run(`read-tree $(s.index)`, dir=dir)      # restore index
    run(`reset -q --soft $(s.head)`, dir=dir) # restore head
end

"""
    Git.transact(f; dir="")

Attempt to execute the function `f`. If this fails, the repository is restored to its
state prior to execution.
"""
function transact(f::Function; dir="")
    state = snapshot(dir=dir)
    try f() catch
        restore(state, dir=dir)
        rethrow()
    end
end

"""
    Git.is_ancestor_of(a, b; dir="")

Determine whether the commit `a` is an ancestor of the commit `b` in the given repository.
"""
function is_ancestor_of(a::AbstractString, b::AbstractString; dir="")
    A = readchomp(`rev-parse $a`, dir=dir)
    readchomp(`merge-base $A $b`, dir=dir) == A
end

const GITHUB_REGEX =
    r"^(?:git@|git://|https://(?:[\w\.\+\-]+@)?)github.com[:/](([^/].+)/(.+?))(?:\.git)?$"i

"""
    Git.set_remote_url(url; remote="origin", dir="")

Add a remote `remote` to the given repository from the URL `url`.
"""
function set_remote_url(url::AbstractString; remote::AbstractString="origin", dir="")
    run(`config remote.$remote.url $url`, dir=dir)
    m = match(GITHUB_REGEX,url)
    m === nothing && return
    push = "git@github.com:$(m.captures[1]).git"
    push != url && run(`config remote.$remote.pushurl $push`, dir=dir)
end

"""
    Git.normalize_url(url)

Normalize the given URL to a valid GitHub repository URL.
"""
function normalize_url(url::AbstractString)
    m = match(GITHUB_REGEX,url)
    m === nothing ? url : "git://github.com/$(m.captures[1]).git"
end

end # module


text:
# This file is a part of Actors.jl, licensed under the MIT License (MIT).

__precompile__()

module Actors

import EasyPkg
EasyPkg.include_all_sources()

function __init__()
end

end # module


text:
module Rock
include("rpc.jl")
include("server.jl")
include("client.jl")
end

max_stars_repo_name: JeremyRubin/Rock


text:
# Koch snowflake
using DrWatson
@quickactivate "NonlinearDynamicsTextbook"
include(srcdir("style.jl"))
using DynamicalSystems, PyPlot, Random

linepoints = SVector{2}.([[0.0; 0.0], [1.0; 0.0]])
flakepoints = SVector{2}.([[0.0; 0.0], [0.5; sqrt(3)/2], [1; 0.0], [0.0; 0.0]])
function koch(points, maxk, α = sqrt(3)/2)
  Q = @SMatrix [0 -1; 1 0]
  for k = 1:maxk
    n = length(points)
    new_points = eltype(points)[]
    for i = 1:n-1
      p1, p2 = points[i], points[i+1]
      v = (p2 - p1) / 3
      q1 = p1 + v
      q2 = p1 + 1.5v + α * Q * v
      q3 = q1 + v
      push!(new_points, p1, q1, q2, q3)
    end
    push!(new_points, points[end])
    points = new_points
  end
  return points
end

# Plot construction
fig, axes = subplots(2,3, figsize = (figx, 2figy))
toprow = axes[1:2:5]
botrow = axes[2:2:6]
kochcolors = COLORS[1:4]

for (o, ax) in enumerate(toprow)
    ax.clear()
    ax.set_aspect("equal")
    k1 = koch(flakepoints, o-1)
    k2 = koch(flakepoints, o)

    ax.plot([a[1] for a in k1], 0.5 .+ [a[2] for a in k1], kochcolors[o], lw = 2.5)
    ax.plot([a[1] for a in k2], 0.5 .+ [a[2] for a in k2], kochcolors[o+1],
    label="Step $(o)", alpha = 0.75, lw = 3)
    # ax.legend(loc = "upper left", framealpha = 1.0, ncol = 2, handlelength = 1)
    ax.text(-0.1, 0.88, "step $(o)", transform = ax.transAxes, size=28)
    # ax.set_ylim(ys...)
    # ax.set_xlim(xs...)
    ax.set_xticklabels([])
    ax.set_yticklabels([])
    ax.axis("off")
end

# Plot versus circle
points1 = [[0.5; 0.865], [2/3, 0.577361]]
points1 = [[0.5; 0.577361], [2/3, 0.865]]
points2 = [[0.537233, 0.801228], [0.555887, 0.768912]]
points2 = [[0.537233, 0.768912], [0.555887,0.801228]]
allpoints = [flakepoints, points1, points2]

for i in 1:3
    botrow[i].clear()
    largekoch = koch(allpoints[i], 6)
    kx, ky = [a[1] for a in largekoch], [a[2] for a in largekoch]
    botrow[i].plot(kx, ky, lw = 1.0)
    botrow[i].set_aspect("equal")
    botrow[i].add_artist(plt.Circle(
        (0.5,0.28900), 0.5, lw=2.0, color="C2", alpha = 0.75, fill=false
    ))
    botrow[i].axis("off")
end

for i in 1:2
    origin = botrow[i]
    zoomin = botrow[i+1]
    zbox = allpoints[i+1]
    axis_zoomin!(zoomin, origin, zbox, zbox, "C$((1,3)[i])"; α = 0.75)
end

fig.tight_layout(pad = 0.1)
fig.subplots_adjust(wspace = 0.01, hspace = 0.01, bottom = 0.01, left = 0.01)
wsave(plotsdir("5", "koch"), fig)

max_stars_repo_name: JuliaDynamics/NonlinearDynamicsTextbook

max_issues_repo_name: JuliaDynamics/NonlinearDynamicsTextbook

max_forks_repo_name: JuliaDynamics/NonlinearDynamicsTextbook


text:
import Dates
import TimeZones

struct WrappedDate
    x::Dates.Date
end
Arrow.ArrowTypes.registertype!(WrappedDate, WrappedDate)

struct WrappedTime
    x::Dates.Time
end
Arrow.ArrowTypes.registertype!(WrappedTime, WrappedTime)

struct WrappedDateTime
    x::Dates.DateTime
end
Arrow.ArrowTypes.registertype!(WrappedDateTime, WrappedDateTime)

struct WrappedZonedDateTime
    x::TimeZones.ZonedDateTime
end
Arrow.ArrowTypes.registertype!(WrappedZonedDateTime, WrappedZonedDateTime)

@testset "Date and time wrappers with missing" begin
    for T in (WrappedDate, WrappedTime, WrappedDateTime, WrappedZonedDateTime)
        if T == WrappedZonedDateTime
            time = T(Dates.now(TimeZones.tz"UTC"))
        else
            time = T(Dates.now())
        end
        table = (; x = [missing, missing, time, missing, time])
        io = IOBuffer()
        Arrow.write(io, table)
        seekstart(io)
        tbl = Arrow.Table(io)
        @test isequal(collect(tbl.x), table.x)
    end
end

@testset "`default(T) isa T`" begin
    for T in (Dates.Date, Dates.Time, Dates.DateTime, TimeZones.ZonedDateTime, Dates.Nanosecond, Dates.Millisecond, Dates.Second, Dates.Day, Dates.Month, Dates.Year)
        @test Arrow.ArrowTypes.default(T) isa T
    end
end


text:
abstract MultivariateFun{T,N}
typealias BivariateFun{T} MultivariateFun{T,2}

export grad, lap, curl

#implements coefficients/values/evaluate
space{T,N}(f::MultivariateFun{T,N})=mapreduce(k->space(f,k),⊗,1:N)
domain{T,N}(f::MultivariateFun{T,N})=mapreduce(k->domain(f,k),*,1:N)

domain(f::MultivariateFun,k::Integer)=domain(space(f,k))

differentiate(u::BivariateFun,i::Integer,j::Integer)=j==0?u:differentiate(differentiate(u,i),i,j-1)
grad(u::BivariateFun)=[differentiate(u,1),differentiate(u,2)]
lap(u::BivariateFun)=differentiate(u,1,2)+differentiate(u,2,2)
Base.div{B<:BivariateFun}(u::Vector{B})=differentiate(u[1],1)+differentiate(u[2],2)
curl{B<:BivariateFun}(u::Vector{B})=differentiate(u[2],1)-differentiate(u[1],2)

Base.chop(f::MultivariateFun)=chop(f,10eps())
Base.eltype{T}(::MultivariateFun{T})=T
Base.eltype{T,N}(::Type{MultivariateFun{T,N}})=T
Base.eltype{MF<:MultivariateFun}(::Type{MF})=eltype(supertype(MF))

include("VectorFun.jl")
include("ProductDomain.jl")
include("TensorSpace.jl")
include("LowRankFun.jl")
include("ProductFun.jl")

arglength(f)=length(Base.uncompressed_ast(f.code.def).args[1])

## Convert between Fun and MultivariateFun
# need to chop extra zeros
Fun(f::ProductFun) =
    Fun(chop!(fromtensor(space(f),coefficients(f)),0),space(f))
Fun(f::ProductFun,sp::TensorSpace) = Fun(ProductFun(f,sp))
Fun(f::LowRankFun) = Fun(ProductFun(f))

Fun(f::MultivariateFun,sp::Space)=Fun(Fun(f),sp)

Fun(f,d1::Domain,d2::Domain)=Fun(f,d1*d2)
Fun{T<:Number,V<:Number}(f,d1::Vector{T},d2::Vector{V})=Fun(f,convert(Domain,d1),convert(Domain,d2))

coefficients(f::BivariateFun,sp::TensorSpace)=coefficients(f,sp[1],sp[2])

points(f::BivariateFun,k...)=points(space(f),size(f,1),size(f,2),k...)

function *(vx::LowRankFun,u0::ProductFun)
    ret=zeros(space(u0))
    for k=1:length(vx.A)
        a,b=vx.A[k],vx.B[k]
        ret+=((b*((a*u0).')).')
    end
    ret
end

*(a::ProductFun,b::LowRankFun)=b*a
*(a::MultivariateFun,b::MultivariateFun)=LowRankFun(a)*ProductFun(b)

for OP in (:+,:-,:*,:/)
    @eval begin
        $OP(f::Fun,g::MultivariateFun)=$OP(ProductFun(f),g)
        $OP(f::MultivariateFun,g::Fun)=$OP(f,ProductFun(g))
    end
end

Base.sum{TS<:TensorSpace}(f::Fun{TS},k::Integer)=sum(ProductFun(f),k)
Base.sum{TS<:TensorSpace}(f::Fun{TS})=sum(ProductFun(f))

## kron

function Base.kron(f::Fun,g::Fun)
    sp=space(f)⊗space(g)
    it=tensorizer(sp)
    N=ncoefficients(f);M=ncoefficients(g)
    cfs=Vector{promote_type(eltype(f),eltype(g))}()
    for (k,j) in it
        # Tensor product is N x M, so if we are outside
        # the (N+M)th diagonal we have no more entries
        if k+j > N+M
            break
        elseif k ≤ N && j ≤ M
            push!(cfs,f.coefficients[k]*g.coefficients[j])
        else
            push!(cfs,0)
        end
    end
    Fun(cfs,sp)
end
Base.kron(f::Fun,g::Number) = kron(f,Fun(g))
Base.kron(f::Number,g::Fun) = kron(Fun(f),g)


text:
module BesselK

  using LinearAlgebra, SpecialFunctions, Polynomials

  export adbesselk, adbesselkxv

  include("besk_ser.jl")   # enhanced direct series. The workhorse for small-ish args.
  include("besk_as.jl")    # asymptotic expansion for large arg.
  include("besk_asv.jl")   # uniform expansion for large order.
  include("besk_temme.jl") # Temme recurrence series for small-ish args. For AD.
  include("besk.jl")       # putting it all together with appropriate branching.

end # module


text:
module Defaults

using FieldMetadata
using FieldMetadata: @default, @redefault, default

export @default_kw, default_kw

macro default_kw(ex)
    typ = get_type(ex)
    quote
        import Defaults.default
        $(FieldMetadata.add_field_funcs(ex, :default))
        $(esc(typ))(;kwargs...) = default_kw($(esc(typ)); kwargs...)
    end
end

macro redefault_kw(ex)
    typ = get_type(ex)
    quote
        import Defaults.default
        $(FieldMetadata.add_field_funcs(ex, :default; update=true))
        $(esc(typ))(;kwargs...) = default_kw($(esc(typ)); kwargs...)
    end
end

get_type(ex) =
    FieldMetadata.firsthead(ex, :struct) do typ_ex
        return FieldMetadata.namify(typ_ex.args[2])
    end

default_kw(::Type{T}; kwargs...) where T = begin
    fnames = fieldnames(T)
    defaults = get_default(T)

    for (key, val) in kwargs
        key in fnames || error("$key is not a field of $T")
        ind = findfirst(n -> n == key, fnames)
        defaults = (defaults[1:ind-1]..., val, defaults[ind+1:end]...)
    end

    T(defaults...)
end

get_default(args...) = default(args...)

end # module


text:
# 64-QAM
"""
---
Calculate a euclidean distance based on input sample and alphabet distribution.
Used for symbolDemappingQAM function
# --- Syntax
	getDist!(out,sR,alphabet, a1...aN)
# --- Input parameters
- out	  : Float64 array of size N
- sR	  : Incoming real sequence
- alphabet : Gray base for decision
- a1...aN : N neigthbours for soft symbol demapper
# --- Output parameters
- []
# ---
# v 1.0 - Robin Gerzaguet.
"""
# 256-QAM
function getDist!(baseComp,sR,alphabet,a1,a2,a3,a4,a5,a6,a7,a8)
	baseComp[1] = abs2.(sR - alphabet[a1]);
	baseComp[2] = abs2.(sR - alphabet[a2]);
	baseComp[3] = abs2.(sR - alphabet[a3]);
	baseComp[4] = abs2.(sR - alphabet[a4]);
	baseComp[5] = abs2.(sR - alphabet[a5]);
	baseComp[6] = abs2.(sR - alphabet[a6]);
	baseComp[7] = abs2.(sR - alphabet[a7]);
	baseComp[8] = abs2.(sR - alphabet[a8]);
	return minimum(baseComp);
end

# 64-QAM
function getDist!(baseComp,sR,alphabet,a1,a2,a3,a4)
	baseComp[1] = abs2.(sR - alphabet[a1]);
	baseComp[2] = abs2.(sR - alphabet[a2]);
	baseComp[3] = abs2.(sR - alphabet[a3]);
	baseComp[4] = abs2.(sR - alphabet[a4]);
	return minimum(baseComp);
end
# 16-QAM
function getDist!(baseComp,sR,alphabet,a1,a2)
	baseComp[1] = abs2.(sR .- alphabet[a1]);
	baseComp[2] = abs2.(sR .- alphabet[a2]);
	return minimum(baseComp);
end

"""
---
Returns the final LLR value based on input distances
# --- Syntax
llr = calcLLR(e0,e1,c)
# --- Input parameters
- e0  : Minimal distance 1
- e1  : Minimal distance 2
- c	  : Channel estimate
# --- Output parameters
- llr : Max likelihood estimate
# ---
# v 1.0 - Robin Gerzaguet.
"""
function calcLLR(e0,e1,c)
	return (e0-e1) * c;
end

"""
---
# --- Description
Returns the log likelihood ratio of the incoming sequence qamSeq based on the channel estimates channelIn. Max log approximation is considered
qamSeq is an input noisy QAM sequence with same size of channel estimate vector
Output is a vector of soft output binary sequence to be fed in a FEC
# --- Syntax
 output	  = symbolDemappingQAM(mcs,qamSeq,channel)
# --- Input parameters
- mcs	  : Constellation size (from 4 to 256) [Int]
- qamSeq  : Complex noisy received sequence (after equalization) [Array{Float64},N]
- channel : Complex channel estimate [Array{Float64},N]
# --- Output parameters
- output  : Soft bits [Array{UInt8},N*log2(mcs)]
"""
function symbolDemappingQAM(mcs,qamSeq,channel)
	# --- Creating output array for allocation
	output = zeros(Float64,Int(length(qamSeq) * log2(mcs)));
	# --- Call bang method
	symbolDemappingQAM!(output,mcs,qamSeq,channel)
	# --- Return array
	return output;
end

"""
---
# --- Description
Returns the log likelihood ratio of the incoming sequence qamSeq based on the channel estimates channelIn
qamSeq is an input noisy QAM sequence with same size of channel estimate vector
Output is populated by  soft output binary sequence to be fed in a FEC
# --- Syntax
 output	  = :symbolDemappingQAM(mcs,qamSeq,channel)
# --- Input parameters
- output  : Soft bits [Array{UInt8},N*log2(mcs)]
- mcs	  : Constellation size (from 4 to 256) [Int]
- qamSeq  : Complex noisy received sequence (after equalization) [Array{Float64},N]
- channel : Complex channel estimate [Array{Float64},N]
# --- Output parameters
- []
# ---
# v 1.0 - Robin Gerzaguet.
"""
function symbolDemappingQAM!(output,mcs,qamSeq,channel)
	# ----------------------------------------------------
	# --- Overall parameters
	# ----------------------------------------------------
	# --- Defining scaling factor
	# Constellation output should have an unitary average power
	# σ^2 = ∑ p(a_i)  ⃒ a_i ⃒ ^2 = 1
	scalingFactor	= sqrt(2/3*(mcs-1));
	# --- Getting size of input
	nbSymb	= length(qamSeq);
	# ----------------------------------------------------
	## --- Symbol demapper
	# ----------------------------------------------------
	if mcs == 4
		# ----------------------------------------------------
		# --- QPSK scheme
		# ----------------------------------------------------
		# Decision is based on euclidean distance with 1/sqrt(2) and -1/sqrt(2).
		# We can do a direct form here.
		output[1:2:end] .= ((abs2.(real(qamSeq).+1/scalingFactor) .- abs2.(real(qamSeq).-1/scalingFactor)).*abs2.(channel))./2;
		output[2:2:end] .= ((abs2.(imag(qamSeq).+1/scalingFactor) .- abs2.(imag(qamSeq).-1/scalingFactor)).*abs2.(channel))./2;
	elseif mcs == 16
		# ----------------------------------------------------
		# --- 16-QAM scheme
		# ----------------------------------------------------
		#maximum log approximation
		# -------------------------%
		#   |      |     |     |
		#  10     11    01    00
		#  -3     -1     1     3
		alphabet  = [-3 -1 1 3]/ scalingFactor;
		# Trying to guess if the first bit  is a 0 or a 1
		# See bitMappingQAM: realSymb = 2bitSeq[4(i-1)+1] + bitSeq[4(i-1)+3]
		# * LSB is third bit for real part
		# * MSB is first bit for real part
		# --> 0 if -3 and 3
		# ---> 1 if -1 and 1
		# Compute here the maxlog approx --> minimum of distance
		# --- We create the container for comparison
		baseComp = zeros(Float64,2);
		# --- Iterative decision
		for iN = 1 : 1 : nbSymb
			c	= abs2(channel[iN]);
			# MSB -> Re
			e_porteur0 = getDist!(baseComp,qamSeq[iN].re,alphabet,3,4)
			e_porteur1 = getDist!(baseComp,qamSeq[iN].re,alphabet,1,2)
			output[4*(iN-1)+1]  = calcLLR(e_porteur0,e_porteur1,c);
			# LSB -> Re
			e_porteur0 = getDist!(baseComp,qamSeq[iN].re,alphabet,1,4)
			e_porteur1 = getDist!(baseComp,qamSeq[iN].re,alphabet,2,3)
			output[4*(iN-1)+3]  =calcLLR(e_porteur0,e_porteur1,c);
			# MSB -> Im
			e_porteur0 = getDist!(baseComp,qamSeq[iN].im,alphabet,3,4)
			e_porteur1 = getDist!(baseComp,qamSeq[iN].im,alphabet,1,2)
			output[4*(iN-1)+2]  =calcLLR(e_porteur0,e_porteur1,c);
			# LSB -> Im
			e_porteur0 = getDist!(baseComp,qamSeq[iN].im,alphabet,1,4)
			e_porteur1 = getDist!(baseComp,qamSeq[iN].im,alphabet,2,3)
			output[4*(iN-1)+4]  =calcLLR(e_porteur0,e_porteur1,c);
		end
	elseif mcs == 64
		# ----------------------------------------------------
		# --- 64-QAM Soft demapper
		# ----------------------------------------------------
		#maximum log approximation
		# ------------------------------------------------#
		#   |      |     |     |     |     |      |      |
		#  100    101   111   110   010   011    001    000
		#   -7     -5    -3    -1    +1    +3     +5     +7
		alphabet  = [-7 -5 -3 -1 1 3 5 7]/ scalingFactor;
		# --- Creating comparison container
		baseComp = zeros(Float64,4);
		# --- Iterative LLR computation
		for iN = 1 : 1 : nbSymb
			c	= abs2(channel[iN]);
			# --- Real part
			# MSB0 - Re
			e_porteur0 = getDist!(baseComp,qamSeq[iN].re,alphabet,5,6,7,8)
			e_porteur1 = getDist!(baseComp,qamSeq[iN].re,alphabet,1,2,3,4)
			output[6*(iN-1)+1]  = calcLLR(e_porteur0,e_porteur1,c);
			# MSB1 - Re
			e_porteur0 = getDist!(baseComp,qamSeq[iN].re,alphabet,1,2,7,8)
			e_porteur1 = getDist!(baseComp,qamSeq[iN].re,alphabet,3,4,5,6)
			output[6*(iN-1)+3]  =calcLLR(e_porteur0,e_porteur1,c);
			# LSB - Re
			e_porteur0 = getDist!(baseComp,qamSeq[iN].re,alphabet,1,4,5,8)
			e_porteur1 = getDist!(baseComp,qamSeq[iN].re,alphabet,2,3,6,7)
			output[6*(iN-1)+5]  =calcLLR(e_porteur0,e_porteur1,c);
			# --- Imag part
			# MSB0 - Im
			e_porteur0 = getDist!(baseComp,qamSeq[iN].im,alphabet,5,6,7,8);
			e_porteur1 = getDist!(baseComp,qamSeq[iN].im,alphabet,1,2,3,4);
			output[6*(iN-1)+2]  = calcLLR(e_porteur0,e_porteur1,c);
			# MSB1 - Im
			e_porteur0 = getDist!(baseComp,qamSeq[iN].im,alphabet,1,2,7,8);
			e_porteur1 = getDist!(baseComp,qamSeq[iN].im,alphabet,3,4,5,6);
			output[6*(iN-1)+4]  = calcLLR(e_porteur0,e_porteur1,c);
			# LSB - Im
			e_porteur0 = getDist!(baseComp,qamSeq[iN].im,alphabet,1,4,5,8);
			e_porteur1 = getDist!(baseComp,qamSeq[iN].im,alphabet,2,3,6,7);
			output[6*(iN-1)+6]  =calcLLR(e_porteur0,e_porteur1,c);
		end
	elseif mcs == 256
		# ----------------------------------------------------
		# --- 256-QAM Soft Demapper
		# ----------------------------------------------------
		#maximum log approximation
		# ------------------------------------------------#
		#	   |       |         |        |        |        |        |        |
		#	0100    0101 .    0111 .   0110 .   0010 .   0011 .   0001 .   0000
		#	 -15     -13 .     -11 .     -9 .     -7 .     -5 .     -3 .     -1
		#	   |       |         |        |        |        |        |        |
		#   1000 .   1001 .   1011 .   1010 .   1110 .   1111 .   1101 .   1100
		#     +1 .     +3 .     +5 .     +7 .     +9 .    +11 .	   +13      +15
		#output    = zeros(Float64,nbSymb*8);
		alphabet  = (-15:2:15)./ scalingFactor;
		# --- Creating comparison container
		baseComp = zeros(Float64,8);
		# --- Iterative LLR computation
		for iN = 1 : 1 : nbSymb
			c	= abs2(channel[iN]);
			# MSB0 - Re
			e_porteur0 = getDist!(baseComp,qamSeq[iN].re,alphabet,1,2,3,4,5,6,7,8);
			e_porteur1 = getDist!(baseComp,qamSeq[iN].re,alphabet,9,10,11,12,13,14,15,16);
			output[8*(iN-1)+1]  = calcLLR(e_porteur0,e_porteur1,c);
			# MSB1 - Re
			e_porteur0 = getDist!(baseComp,qamSeq[iN].re,alphabet,5,6,7,8,9,10,11,12);
			e_porteur1 = getDist!(baseComp,qamSeq[iN].re,alphabet,1,2,3,4,13,14,15,16);
			output[8*(iN-1)+3]  = calcLLR(e_porteur0,e_porteur1,c);
			# MSB2 - Re
			e_porteur0 = getDist!(baseComp,qamSeq[iN].re,alphabet,1,2,7,8,9,10,15,16);
			e_porteur1 = getDist!(baseComp,qamSeq[iN].re,alphabet,3,4,5,6,11,12,13,14);
			output[8*(iN-1)+5]  = calcLLR(e_porteur0,e_porteur1,c);
			# LSB - Re
			e_porteur0 = getDist!(baseComp,qamSeq[iN].re,alphabet,1,4,5,8,9,12,13,16);
			e_porteur1 = getDist!(baseComp,qamSeq[iN].re,alphabet,2,3,6,7,10,11,14,15);
			output[8*(iN-1)+7]  = calcLLR(e_porteur0,e_porteur1,c);
			# --- Imag part config
			# MSB0 - Im
			e_porteur0 = getDist!(baseComp,qamSeq[iN].im,alphabet,1,2,3,4,5,6,7,8);
			e_porteur1 = getDist!(baseComp,qamSeq[iN].im,alphabet,9,10,11,12,13,14,15,16);
			output[8*(iN-1)+2]  = calcLLR(e_porteur0,e_porteur1,c);
			# MSB1 - Im
			e_porteur0 = getDist!(baseComp,qamSeq[iN].im,alphabet,5,6,7,8,9,10,11,12);
			e_porteur1 = getDist!(baseComp,qamSeq[iN].im,alphabet,1,2,3,4,13,14,15,16);
			output[8*(iN-1)+4]  = calcLLR(e_porteur0,e_porteur1,c);
			# MSB2 - Im
			e_porteur0 = getDist!(baseComp,qamSeq[iN].im,alphabet,1,2,7,8,9,10,15,16);
			e_porteur1 = getDist!(baseComp,qamSeq[iN].im,alphabet,3,4,5,6,11,12,13,14);
			output[8*(iN-1)+6]  = calcLLR(e_porteur0,e_porteur1,c);
			# LSB - Im
			e_porteur0 = getDist!(baseComp,qamSeq[iN].im,alphabet,1,4,5,8,9,12,13,16);
			e_porteur1 = getDist!(baseComp,qamSeq[iN].im,alphabet,2,3,6,7,10,11,14,15);
			output[8*(iN-1)+8]  = calcLLR(e_porteur0,e_porteur1,c);
		end
	end
	return output
end

# ----------------------------------------------------
# --- Convert LLR to Unsigned Char Array (soft decoder in C)
# ----------------------------------------------------
"""
---
Convert a LLR floating array to a UInt LLR
Some method (espcially C functions as in libfec) expect LLR to be UInt value from  0 (Likely a 0) to 255 (Likely a 1).
symbolDemappingQAM create a Float64 array from -infty (likely a 0) to nfty (Likely a 1). This function take a Float64 LLR estimate and output a UInt8 vector
# --- Syntax
	llrToUInt!(llrUInt,llr)
# --- Input parameters
- llrUInt : Output LLR in UInt [Array{UInt},N]
- llr	  : Input LLR [Array{Float64},N]
# --- Output parameters
- []
# ---
# v 1.0 - Robin Gerzaguet.
"""
function llrToUInt!(llrUInt,llr)
    # 0 ---> Likely a 0
    # 255 ---> Likely a 1
    # We have the llRin: >0 1 <0 -1, higher the better reliability.
    # We assume that 1  --> certainly a 1 (i.e 255)
    #                -1 --> certainly a 0 (i.e 0)
    # We need to convert LLR into UInt8: LLR is converted between 0 and 1 and map to 256 values
    # 0 is a 0 with high probability
    # 1 is a 1 with high probability
    # --- Minimal and maximal bounds (-1 +1)
    # Then, Switch to [0;2]
    # Then Convert to [0,1]
	# Finally extend to 255
	llrUInt .= UInt8.( round.(255* (1 .+ map(x->min(x,1),map(x->max(x,-1),llr)))/2));
end

"""
---
Convert a LLR floating array to a UInt LLR
Some method (espcially C functions as in libfec) expect LLR to be UInt value from  0 (Likely a 0) to 255 (Likely a 1).
symbolDemappingQAM create a Float64 array from -infty to infty. This function take a Float64 LLR estimate and output a UInt8 vector
# --- Syntax
	llrUInt = llrToUInt(llr)
# --- Input parameters
- llr	  : Input LLR [Array{Float64},N]
# --- Output parameters
- llrUInt : Output LLR in UInt [Array{UInt},N]
# ---
# v 1.0 - Robin Gerzaguet.
"""
function llrToUInt(llr)
	# --- Create output vector
	llrUInt = zeros(UInt8,length(llr));
	# --- Call bang method
	llrToUInt!(llrUInt,llr);
	# --- Return LLR
    return llrUInt;
end

# ----------------------------------------------------
# --- Hard decoding from Soft LLRs
# ----------------------------------------------------
"""
---
Returns hard binary value from soft LLR estimate (no FEC decoder, only hard decision here!)
# --- Syntax
	llrToHardBits(hardD,llr)
# --- Input parameters
- hardD	  : Hard binary decision [Array{UInt},N]
- llr	  : Input LLR (Float64 or UInt8 array) [Union{[Array{UInt}],[Array{Float64}]}]
# --- Output parameters
- []
# ---
# v 1.0 - Robin Gerzaguet.
"""
function llrToHardBits!(hardD,llr::Array{Float64})
	# --- Direct LLR are given, sign is bit information
	hardD   .= map(x-> (x<0) ? 0 : 1,llr);
	return hardD
end

"""
---
Returns hard binary value from soft LLR estimate (no FEC decoder, only hard decision here!)
# --- Syntax
	llrToHardBits(hardD,llr)
# --- Input parameters
- hardD	  : Hard binary decision [Array{UInt},N]
- llr	  : Input LLR (Float64 or UInt8 array) [Union{[Array{UInt}],[Array{Float64}]}]
# --- Output parameters
- []
# ---
# v 1.0 - Robin Gerzaguet.
"""
function llrToHardBits!(hardD,llr::Array{UInt8})
	# --- UInt LLR (input of soft decoder): rule is 127
	hardD   .= map(x-> (x<0x7f) ? 0 : 1,llr);
	return hardD
end

"""
---
Returns hard binary value from soft LLR estimate (no FEC decoder, only hard decision here!)
# --- Syntax
	hardD = llrToHardBits(llr)
# --- Input parameters
- llr	  : Input LLR (Float64 or UInt8 array) [Union{[Array{UInt}],[Array{Float64}]}]
# --- Output parameters
- hardD	  : Hard binary decision [Array{UInt},N]
# ---
# v 1.0 - Robin Gerzaguet.
"""
function llrToHardBits(llr::Array{Float64})
	# --- Direct LLR are given, sign is bit information
	hardD   = map(x-> (x<0) ? 0 : 1,llr);
	return hardD
end

"""
---
Returns hard binary value from soft LLR estimate (no FEC decoder, only hard decision here!)
# --- Syntax
	hardD = llrToHardBits(llr)
# --- Input parameters
- llr	  : Input LLR (Float64 or UInt8 array) [Union{[Array{UInt}],[Array{Float64}]}]
# --- Output parameters
- hardD	  : Hard binary decision [Array{UInt},N]
# ---
# v 1.0 - Robin Gerzaguet.
"""
function llrToHardBits(llr::Array{UInt8})
	# --- UInt LLR (input of soft decoder): rule is 127
	hardD   = map(x-> (x<0x7f) ? 0 : 1,llr);
	return hardD
end


text:
function update_omega!(state::State, data::Data, prior::Prior)
  a, b = prior.a_omega, state.tau
  akernel = zero.(state.omega)
  bkernel = zero.(state.omega)

  for i in samplenames
    y = ref_yfinite(data, i)
    lam = ref_lambda(state, i)
    zeta = ref_zeta(state, i)
    v = ref_v(state, i)
    for n in eachindex(y)
      k = lam[n]
      akernel[k] += 1
      bkernel[k] += v[n] * (y[n] - state.mu[k] - state.psi[k] * zeta[n])^2
    end
  end

  anew = a .+ akernel / 2
  bnew = b .+ bkernel / 2

  for k in 1:prior.K
    state.omega[k] = rand(InverseGamma(anew[k], bnew[k]))
  end

  update_sigma!(state)
  update_phi!(state)
end

max_stars_repo_name: luiarthur/CytofDensityEstimation

max_issues_repo_name: luiarthur/CytofDensityEstimation

max_forks_repo_name: luiarthur/CytofDensityEstimation


text:
module ReactionMechanismSimulator
    using PyCall
    push!(PyVector(pyimport("sys")["path"]), "")
    const Chem = PyNULL()
    const Desc = PyNULL()
    const molecule = PyNULL()
    const pydot = PyNULL()
    const chemkin =PyNULL()
    const species = PyNULL()
    const reaction = PyNULL()
    const nasa = PyNULL()
    const wilhoit = PyNULL()
    const arrhenius = PyNULL()
    const falloff = PyNULL()
    const chebyshev = PyNULL()
    const solvation = PyNULL()
    const yaml = PyNULL()
    const os = PyNULL()
    function __init__()
        copy!(Chem,pyimport_conda("rdkit.Chem","rdkit","rdkit"))
        copy!(Desc,pyimport_conda("rdkit.Chem.Descriptors","rdkit","rdkit"))
        copy!(molecule,pyimport_conda("rmgpy.molecule","rmg","rmg"))
        copy!(pydot,pyimport_conda("pydot","pydot"))
        copy!(chemkin,pyimport_conda("rmgpy.chemkin","rmg","rmg"))
        copy!(species,pyimport_conda("rmgpy.species","rmg","rmg"))
        copy!(reaction,pyimport_conda("rmgpy.reaction","rmg","rmg"))
        copy!(nasa,pyimport_conda("rmgpy.thermo.nasa","rmg","rmg"))
        copy!(wilhoit,pyimport_conda("rmgpy.thermo.wilhoit","rmg","rmg"))
        copy!(arrhenius,pyimport_conda("rmgpy.kinetics.arrhenius","rmg","rmg"))
        copy!(falloff,pyimport_conda("rmgpy.kinetics.falloff","rmg","rmg"))
        copy!(chebyshev,pyimport_conda("rmgpy.kinetics.chebyshev","rmg","rmg"))
        copy!(solvation,pyimport_conda("rmgpy.data.solvation","rmg","rmg"))
        copy!(yaml,pyimport_conda("yaml","yaml"))
        copy!(os,pyimport_conda("os","os"))
    end
    include("Constants.jl")
    include("Tools.jl")
    include("Calculators/RateUncertainty.jl")
    include("Calculators/ThermoUncertainty.jl")
    include("Calculators/Thermo.jl")
    include("Calculators/Diffusion.jl")
    include("Calculators/Rate.jl")
    include("Calculators/Viscosity.jl")
    include("Calculators/Thermovec.jl")
    include("Calculators/Ratevec.jl")
    include("Calculators/kLAkH.jl")
    include("Species.jl")
    include("Solvent.jl")
    include("Reaction.jl")
    include("Phase.jl")
    include("PhaseState.jl")
    include("Interface.jl")
    include("Domain.jl")
    include("yml.jl")
    include("Parse.jl")
    include("Reactor.jl")
    include("Simulation.jl")
    include("TransitorySensitivities.jl")
    include("AutomaticMechanismAnalysis.jl")
    include("EdgeAnalysis.jl")
    include("Debugging.jl")
    include("Plotting.jl")
    include("fluxdiagrams.jl")
end


text:
using DiffEqFlux, OrdinaryDiffEq, Test

function f(u, p, t)
    p .* u
end
rc=62
ps = repeat([-0.001], rc)
tspan = (7.0, 84.0)
u0 = 3.4 .+ ones(rc)
t = collect(range(minimum(tspan), stop=maximum(tspan), length=157))
prob = ODEProblem(f, u0, tspan, ps)
data = Array(solve(prob, Tsit5(), saveat = t, abstol=1e-12, reltol=1e-12))
ptest = ones(rc)

u′,u = collocate_data(data,t,SigmoidKernel())
@test sum(abs2,u - data) < 1e-8

function loss(p)
    cost = zero(first(p))
    for i in 1:size(u′,2)
      _du = f(@view(u[:,i]),p,t[i])
      u′i = @view u′[:,i]
      cost += sum(abs2,u′i .- _du)
    end
    sqrt(cost)
end
@test loss(ptest) ≈ 418.3400017500223


text:
import accrue.cryptoerase.runtime.Condition;

class H {
    public C c;
    public H(C c) {
	this.c = c;
    }
}

class C {
    static final Condition c = new Condition();
    int{L /c H} x = 0;

    void put(int v) {
	this.x = v;
    }

    int get() {
	return this.x;
    }

    public static void main(String[] args) {
	C c = new C();
	H h = new H(c);
	c.put(42);
	c.put(h.c.get());
    }
}

max_stars_repo_name: HarvardPL/cryptoerase

max_issues_repo_name: HarvardPL/cryptoerase

max_forks_repo_name: HarvardPL/cryptoerase


text:
@testset "ColSpec" begin
  vecnames = [:a, :b, :c, :d, :e, :f]
  tupnames = (:a, :b, :c, :d, :e, :f)

  # vector of symbols
  colspec = [:a, :c, :e]
  snames = TableTransforms.choose(colspec, vecnames)
  @test snames == [:a, :c, :e]
  snames = TableTransforms.choose(colspec, tupnames)
  @test snames == [:a, :c, :e]

  # tuple of symbols
  colspec = (:a, :c, :e)
  snames = TableTransforms.choose(colspec, vecnames)
  @test snames == [:a, :c, :e]
  snames = TableTransforms.choose(colspec, tupnames)
  @test snames == [:a, :c, :e]

  # vector of strings
  colspec = ["a", "c", "e"]
  snames = TableTransforms.choose(colspec, vecnames)
  @test snames == [:a, :c, :e]
  snames = TableTransforms.choose(colspec, tupnames)
  @test snames == [:a, :c, :e]

  # tuple of strings
  colspec = ("a", "c", "e")
  snames = TableTransforms.choose(colspec, vecnames)
  @test snames == [:a, :c, :e]
  snames = TableTransforms.choose(colspec, tupnames)
  @test snames == [:a, :c, :e]

  # vector of integers
  colspec = [1, 3, 5]
  snames = TableTransforms.choose(colspec, vecnames)
  @test snames == [:a, :c, :e]
  snames = TableTransforms.choose(colspec, tupnames)
  @test snames == [:a, :c, :e]

  # tuple of integers
  colspec = (1, 3, 5)
  snames = TableTransforms.choose(colspec, vecnames)
  @test snames == [:a, :c, :e]
  snames = TableTransforms.choose(colspec, tupnames)
  @test snames == [:a, :c, :e]

  # regex
  colspec = r"[ace]"
  snames = TableTransforms.choose(colspec, vecnames)
  @test snames == [:a, :c, :e]
  snames = TableTransforms.choose(colspec, tupnames)
  @test snames == [:a, :c, :e]

  # colon
  snames = TableTransforms.choose(:, vecnames)
  @test snames == [:a, :b, :c, :d, :e, :f]
  snames = TableTransforms.choose(:, tupnames)
  @test snames == [:a, :b, :c, :d, :e, :f]

  # throws
  @test_throws AssertionError TableTransforms.choose(r"x", vecnames)
  @test_throws AssertionError TableTransforms.choose(r"x", tupnames)
  @test_throws AssertionError TableTransforms.choose(String[], vecnames)
  @test_throws AssertionError TableTransforms.choose(String[], tupnames)
  @test_throws AssertionError TableTransforms.choose(Symbol[], vecnames)
  @test_throws AssertionError TableTransforms.choose(Symbol[], tupnames)

  # type stability
  @inferred TableTransforms.choose([:a, :b], vecnames)
  @inferred TableTransforms.choose([:a, :b], tupnames)
  @inferred TableTransforms.choose((:a, :b), vecnames)
  @inferred TableTransforms.choose((:a, :b), tupnames)
  @inferred TableTransforms.choose(["a", "b"], vecnames)
  @inferred TableTransforms.choose(["a", "b"], tupnames)
  @inferred TableTransforms.choose(("a", "b"), vecnames)
  @inferred TableTransforms.choose(("a", "b"), tupnames)
  @inferred TableTransforms.choose([1, 2], vecnames)
  @inferred TableTransforms.choose([1, 2], tupnames)
  @inferred TableTransforms.choose((1, 2), vecnames)
  @inferred TableTransforms.choose((1, 2), tupnames)
  @inferred TableTransforms.choose(r"[ab]", vecnames)
  @inferred TableTransforms.choose(r"[ab]", tupnames)
  @inferred TableTransforms.choose(:, vecnames)
  @inferred TableTransforms.choose(:, tupnames)
end


text:
using Media, Lazy, Hiccup

import Media: render, @render
import Hiccup: div

import Juno: Inline, Clipboard, Editor, Console, PlotPane, Model, Tree, LazyTree,
             SubTree, Copyable, Link, link, Row, Table, fade, icon, interleave, dims,
             undefs, errtrace, errmsg, view, UNDEF, limit

Media.@defpool Editor
Media.@defpool Console

setdisplay(Editor(), Any, Console())
setdisplay(Console(), Any, Console())

render(::Clipboard, x) =
  sprint(io -> show(IOContext(io, :limit => true), MIME"text/plain"(), x))

render(::Editor, x) =
  render(Inline(), Copyable(x))

render(::Console, x) =
  Atom.msg("result", render(Inline(), Copyable(x)))

render(::PlotPane, x) =
  Atom.msg("plot", render(Inline(), x))

@render Inline l::Row begin
  span([render(Inline(), x) for x in l.xs])
end

include("showdisplay.jl")
include("base.jl")
include("utils.jl")
include("view.jl")
include("lazy.jl")
include("errors.jl")
include("frontend.jl")


text:
function hasfieldcount(@nospecialize(dt))
    try
        fieldcount(dt)
    catch
        return false
    end
    return true
end

if VERSION <= v"1.6"
    allocatedinline(@nospecialize(T)) = T.isinlinealloc
else
    import Base: allocatedinline
end


text:
# Weave readme
using Pkg
cd("c:/git/JDF/readme-build")
Pkg.activate("c:/git/JDF/readme-build")
Pkg.update()
upcheck()

using Weave

weave("../README.jmd", out_path="c:/git/JDF/", doctype="github")

if false
    # debug
    tangle("../README.jmd")
end

max_stars_repo_name: xiaodaigh/JDF

max_issues_repo_name: xiaodaigh/JDF

max_forks_repo_name: xiaodaigh/JDF


text:
import Random: rand, rand!

abstract type VSLDiscreteDistribution{T} <: VSLDistribution{T} end

macro vsl_distribution_discrete(name, methods, rtype, properties...)
    t = :(Union{})
    method_symbols = [Symbol(string("VSL_RNG_METHOD_", method)) for method in methods.args]
    for method in method_symbols
        push!(t.args, method)
    end
    code = quote
        mutable struct $name{T<:$t} <: VSLDiscreteDistribution{$rtype}
            brng::BasicRandomNumberGenerator
            ii::Int
            tmp::Vector{$rtype}
        end
    end
    fields = code.args[2].args[end].args
    for property in properties
        push!(fields, property)
    end
    push!(fields, :(method::Type{T}))

    push!(code.args, :(
        $name(brng::BasicRandomNumberGenerator, $(properties...), method::Type{T}) where {T<:$t} =
        $name(brng, BUFFER_LENGTH, zeros($rtype, BUFFER_LENGTH), $([property.args[1] for property in properties]...), method)
    ))

    push!(code.args, :(
        $name(brng::BasicRandomNumberGenerator, $(properties...)) =
        $name(brng, $([property.args[1] for property in properties]...), $(method_symbols[1]))
    ))

    esc(code)
end

macro register_rand_functions_discrete(typename, name, methods, method_constants, rtype, types, arguments)
    code = quote end
    method_symbols = [Symbol(string("VSL_RNG_METHOD_", method)) for method in methods.args]
    function_name = "viRng$(name)"
    for (method, constant) in zip(method_symbols, method_constants.args)
        push!(
            code.args,
            :(function rand(d::$typename{$method}, ::Type{$rtype}=$rtype)
                d.ii += 1
                if d.ii > BUFFER_LENGTH
                    ccall(($function_name, libmkl), Cint, (Int, Ptr{Cvoid}, Int, Ptr{$rtype}, $(types.args...)),
                        $constant, d.brng.stream_state[1], BUFFER_LENGTH, d.tmp, $(arguments.args...))
                    d.ii = 1
                end
                d.tmp[d.ii]
            end)
        )
        push!(
            code.args,
            :(function rand!(d::$typename{$method}, A::Array{$rtype}, ::Type{$rtype}=$rtype)
                n = length(A)
                t = BUFFER_LENGTH - d.ii
                if n <= t
                    copyto!(A, 1, d.tmp, d.ii + 1, n)
                    d.ii = d.ii + n
                    return A
                end
                copyto!(A, 1, d.tmp, d.ii + 1, t)
                d.ii = BUFFER_LENGTH
                ccall(($function_name, libmkl), Cint, (Int, Ptr{Cvoid}, Int, Ptr{$rtype}, $(types.args...)),
                    $constant, d.brng.stream_state[1], n - t, view(A, t+1:n), $(arguments.args...))
                A
            end)
        )
        push!(
            code.args,
            :(rand(d::$typename{$method}, dims::Dims) = rand!(d, Array{$rtype}(dims)))
        )
    end
    esc(code)
end

@vsl_distribution_discrete(
    UniformDiscrete,
    [STD],
    Cint,
    a::Cint,
    b::Cint
)

@register_rand_functions_discrete(
    UniformDiscrete,
    Uniform,
    [STD],
    [0x00000000],
    Cint,
    [Cint, Cint],
    [d.a, d.b]
)

@vsl_distribution_discrete(
    UniformBits,
    [STD],
    Cuint,
)

@register_rand_functions_discrete(
    UniformBits,
    UniformBits,
    [STD],
    [0x00000000],
    Cuint,
    [],
    []
)

@vsl_distribution_discrete(
    UniformBits32,
    [STD],
    Cuint,
)

@register_rand_functions_discrete(
    UniformBits32,
    UniformBits32,
    [STD],
    [0x00000000],
    Cuint,
    [],
    []
)

@vsl_distribution_discrete(
    UniformBits64,
    [STD],
    UInt64,
)

@register_rand_functions_discrete(
    UniformBits64,
    UniformBits64,
    [STD],
    [0x00000000],
    UInt64,
    [],
    []
)

@vsl_distribution_discrete(
    Bernoulli,
    [ICDF],
    Cint,
    p::Cdouble
)

@register_rand_functions_discrete(
    Bernoulli,
    Bernoulli,
    [ICDF],
    [0x00000000],
    Cint,
    [Cdouble],
    [d.p]
)

@vsl_distribution_discrete(
    Geometric,
    [ICDF],
    Cint,
    p::Cdouble
)

@register_rand_functions_discrete(
    Geometric,
    Geometric,
    [ICDF],
    [0x00000000],
    Cint,
    [Cdouble],
    [d.p]
)

@vsl_distribution_discrete(
    Binomial,
    [BTPE],
    Cint,
    ntrial::Cint,
    p::Cdouble
)

@register_rand_functions_discrete(
    Binomial,
    Binomial,
    [BTPE],
    [0x00000000],
    Cint,
    [Cint, Cdouble],
    [d.ntrial, d.p]
)

@vsl_distribution_discrete(
    Hypergeometric,
    [H2PE],
    Cint,
    l::Cint,
    s::Cint,
    m::Cint
)

@register_rand_functions_discrete(
    Hypergeometric,
    Hypergeometric,
    [H2PE],
    [0x00000000],
    Cint,
    [Cint, Cint, Cint],
    [d.l, d.s, d.m]
)

@vsl_distribution_discrete(
    Poisson,
    [PTPE, POISNORM],
    Cint,
    λ::Cdouble
)

@register_rand_functions_discrete(
    Poisson,
    Poisson,
    [PTPE, POISNORM],
    [0x00000000, 0x00000001],
    Cint,
    [Cdouble],
    [d.λ]
)

@vsl_distribution_discrete(
    PoissonV,
    [POISNORM],
    Cint,
    λ::Vector{Cdouble}
)

@register_rand_functions_discrete(
    PoissonV,
    PoissonV,
    [POISNORM],
    [0x00000000],
    Cint,
    [Ptr{Cdouble}],
    [d.λ]
)

@vsl_distribution_discrete(
    NegBinomial,
    [NBAR],
    Cint,
    a::Cdouble,
    p::Cdouble
)

@register_rand_functions_discrete(
    NegBinomial,
    NegBinomial,
    [NBAR],
    [0x00000000],
    Cint,
    [Cdouble, Cdouble],
    [d.a, d.p]
)


text:
# Problem 2 - Project Euler
# http://projecteuler.net/index.php?section=problems&id=2
function fibevensum(a, b, sum, xmax)
    if a >= xmax
        sum
    elseif a % 2 == 0
        fibevensum(b, a + b, sum + a, xmax)
    else
        fibevensum(b, a + b, sum, xmax)
    end
end

println(fibevensum(1,2,0, 4000000))

max_stars_repo_name: emergent/ProjectEuler

max_issues_repo_name: emergent/ProjectEuler

max_forks_repo_name: emergent/ProjectEuler


text:
# The following variables need to be defined before `include`-ing this file:
# `testset_name_suffix`
# `n_values`
# `k_values`
# `m_values`
function test_complex(::Type{TE}, m_values, k_values, n_values, testset_name_suffix) where {TE}
  T = Complex{TE}
  @time @testset "Matrix Multiply $T $(testset_name_suffix)" begin
    for n ∈ n_values
      for k ∈ k_values
        for m ∈ m_values
          A = rand(T, m, k);
          B = rand(T, k, n);
          b = rand(T, k);

          Are = real.(A);
          Bre = real.(B);
          bre = real.(b);

          A′ = permutedims(A)'
          B′ = permutedims(B)'
          AB = A * B;
          A′B = A′*B;
          AB′ = A*B′;
          A′B′= A′*B′;
          Ab = A*b;
          A′b = A′*b;

          AreB = Are*B;
          ABre = A*Bre;
          Areb = Are*b;
          Abre = A*bre;

          @info "" T n k m
          @test @time(Octavian.matmul(A, B)) ≈ AB
          @test @time(Octavian.matmul(A, Bre)) ≈ ABre
          @test @time(Octavian.matmul(Are, B)) ≈ AreB
          @test @time(Octavian.matmul(A′, B)) ≈ A′B
          @test @time(Octavian.matmul(A, B′)) ≈ AB′
          @test @time(Octavian.matmul(A′, B′)) ≈ A′B′

          @test @time(Octavian.matmul(A, b)) ≈ Ab
          @test transpose(@time(Octavian.matmul(transpose(b), transpose(A)))) ≈ Ab
          @test @time(Octavian.matmul(A, bre)) ≈ Abre
          @test @time(Octavian.matmul(Are, b)) ≈ Areb
          @test @time(Octavian.matmul(A′, b)) ≈ A′b
          @test transpose(@time(Octavian.matmul(transpose(b), transpose(A′)))) ≈ A′b

          @test @time(Octavian.matmul_serial(A, B)) ≈ AB
          @test @time(Octavian.matmul_serial(A, Bre)) ≈ ABre
          @test @time(Octavian.matmul_serial(Are, B)) ≈ AreB
          @test @time(Octavian.matmul_serial(A′, B)) ≈ A′B
          @test @time(Octavian.matmul_serial(A, B′)) ≈ AB′
          @test @time(Octavian.matmul_serial(A′, B′)) ≈ A′B′

          @test @time(Octavian.matmul_serial(A, b)) ≈ Ab
          @test transpose(@time(Octavian.matmul_serial(transpose(b), transpose(A)))) ≈ Ab
          @test @time(Octavian.matmul_serial(A, bre)) ≈ Abre
          @test @time(Octavian.matmul_serial(Are, b)) ≈ Areb
          @test @time(Octavian.matmul_serial(A′, b)) ≈ A′b
          @test transpose(@time(Octavian.matmul_serial(transpose(b), transpose(A′)))) ≈ A′b

          C = Matrix{T}(undef, n, m)'
          @test @time(Octavian.matmul!(C, A, B)) ≈ AB

          C1 = rand(T, m, n)
          C2 = copy(C1)
          α, β = T(1 - 2im), T(3 + 4im)
          @test @time(Octavian.matmul!(C1, A, B, α, β)) ≈ Octavian.matmul!(C2, A, B, α, β)
        end
      end
    end
  end
end

function matmul_pack_ab!(C, A, B)
  M, N = size(C); K = size(B,1)
  zc, za, zb = Octavian.zstridedpointer.((C,A,B))
  nspawn = min(Threads.nthreads(), Octavian.num_cores())
  GC.@preserve C A B begin
    if nspawn > 1
      threads, torelease = Octavian.PolyesterWeave.__request_threads((nspawn-1)%UInt32, Octavian.PolyesterWeave.worker_pointer(), nothing)
      @assert threads.i < Threads.nthreads()
      Octavian.matmul_pack_A_and_B!(
        zc, za, zb, Octavian.StaticInt{1}(), Octavian.StaticInt{0}(), M, K, N, threads,
        Octavian.W₁Default(), Octavian.W₂Default(), Octavian.R₁Default(), Octavian.R₂Default()
      )
      Octavian.PolyesterWeave.free_threads!(torelease)
    else
      Octavian.matmul_st_pack_A_and_B!(
        zc, za, zb, Octavian.StaticInt{1}(), Octavian.StaticInt{0}(), M, K, N,
        Octavian.W₁Default(), Octavian.W₂Default(), Octavian.R₁Default(), Octavian.R₂Default(), 0
      )
    end
  end
  C
end

function test_real(::Type{T}, m_values, k_values, n_values, testset_name_suffix) where {T}
  @time @testset "Matrix Multiply $T $(testset_name_suffix)" begin
    for n ∈ n_values
      for k ∈ k_values
        for m ∈ m_values
          A = rand(T, m, k)
          B = rand(T, k, n)
          b = rand(T, k)
          A′ = permutedims(A)'
          B′ = permutedims(B)'
          AB = A * B;
          Ab = A * b;
          @info "" T n k m
          @test @time(Octavian.matmul(A, B)) ≈ AB
          @test @time(Octavian.matmul(A′, B)) ≈ AB
          @test @time(Octavian.matmul(A, B′)) ≈ AB
          @test @time(Octavian.matmul(A′, B′)) ≈ AB
          @test @time(Octavian.matmul_serial(A, B)) ≈ AB
          @test @time(Octavian.matmul_serial(A′, B)) ≈ AB
          @test @time(Octavian.matmul_serial(A, B′)) ≈ AB
          @test @time(Octavian.matmul_serial(A′, B′)) ≈ AB
          @test @time(Octavian.matmul(A, b)) ≈ Ab
          @test @time(Octavian.matmul(A′, b)) ≈ Ab
          @test @time(Octavian.matmul(b', A'))' ≈ Ab
          @test @time(Octavian.matmul(b', A′'))' ≈ Ab
          @test @time(Octavian.matmul_serial(A, b)) ≈ Ab
          @test @time(Octavian.matmul_serial(A′, b)) ≈ Ab
          @test @time(Octavian.matmul_serial(b', A'))' ≈ Ab
          @test @time(Octavian.matmul_serial(b', A′'))' ≈ Ab
        end
      end
    end
    m = k = n = max(8Octavian.OCTAVIAN_NUM_TASKS[], 400)
    A = rand(T, m, k);
    B = rand(T, k, n);
    A′ = permutedims(A)';
    B′ = permutedims(B)';
    AB = A * B;
    @test matmul_pack_ab!(similar(AB), A, B) ≈ AB
    @test matmul_pack_ab!(similar(AB), A, B′) ≈ AB
    @test matmul_pack_ab!(similar(AB), A′, B) ≈ AB
    @test matmul_pack_ab!(similar(AB), A′, B′) ≈ AB
  end
end

@time @testset "zero-sized-matrices" begin
    @test Octavian.matmul_serial(randn(0,0), randn(0,0)) == zeros(0, 0)
    @test Octavian.matmul_serial(randn(2,3), randn(3,0)) == zeros(2, 0)
    @test Octavian.matmul_serial(randn(2,0), randn(0,2)) == zeros(2, 2)
    @test Octavian.matmul_serial!(ones(2,2),randn(2,0), randn(0,2), 1.0, 2.0) == ones(2, 2) .* 2
    @test Octavian.matmul(randn(0,0), randn(0,0)) == zeros(0, 0)
    @test Octavian.matmul(randn(2,3), randn(3,0)) == zeros(2, 0)
    @test Octavian.matmul(randn(2,0), randn(0,2)) == zeros(2, 2)
    @test Octavian.matmul!(ones(2,2),randn(2,0), randn(0,2), 1.0, 2.0) == ones(2, 2) .* 2
end


text:
module ParametricCVforSFS

export parametric_SFS_SI_TVS, parametric_SFS_SI_CV, parametric_SFS_CI_TVS, parametric_SFS_CI_CV

include("../data.jl")
using .DataUtil
include("sfs.jl")
using .SFS
include("selective_inference.jl")
using .SIforSFS
include("../quadratic.jl")
using .Quadratic
include("../intervals.jl")
using .Intersection
include("../selective_p.jl")
using .MultiTruncatedDistributions
include("../ci.jl")
using .ConfidenceInterval
using IntervalSets
using Distributions

const ε = 1e-4

function validation_error(𝐗_train, 𝐗_val, 𝐚_train, 𝐚_val, 𝐛_train, 𝐛_val)
    𝐗₊ = (𝐗_train' * 𝐗_train) \ 𝐗_train' # pseudo inverse matrix of 𝐗_train
    𝐚 = 𝐚_val .- 𝐗_val * 𝐗₊ * 𝐚_train
    𝐛 = 𝐛_val .- 𝐗_val * 𝐗₊ * 𝐛_train
    return quadratic(𝐛' * 𝐛, 2 * 𝐚' * 𝐛, 𝐚' * 𝐚) # E(z) = 𝐚ᵀ𝐚 + 2𝐚ᵀ𝐛z + 𝐛ᵀ𝐛z²
end

function compute_val_error_path(k, 𝐗_train, 𝐗_val, 𝐚_train, 𝐚_val, 𝐛_train, 𝐛_val, z_min, z_max)
    𝐳_k, 𝑨_train_k, _ = compute_solution_path(𝐗_train, 𝐚_train, 𝐛_train, k, z_min, z_max)
    𝐄_k = (𝐴 -> validation_error((@view 𝐗_train[:, 𝐴]), (@view 𝐗_val[:, 𝐴]), 𝐚_train, 𝐚_val, 𝐛_train, 𝐛_val)).(𝑨_train_k)
    return 𝐳_k, 𝐄_k
end

function compute_cv_path(k, 𝐗, 𝐚, 𝐛, 𝑘, z_min, z_max)
    𝑘_fold = KFold(length(𝐚), 𝑘) # 𝑘-fold Cross Validation
    paths = (i -> begin
        𝐗_train, 𝐗_val = train_val_split(𝐗, 𝑘_fold.parts[i], 𝑘_fold.parts[i + 1] - 1)
        𝐚_train, 𝐚_val = train_val_split(𝐚, 𝑘_fold.parts[i], 𝑘_fold.parts[i + 1] - 1)
        𝐛_train, 𝐛_val = train_val_split(𝐛, 𝑘_fold.parts[i], 𝑘_fold.parts[i + 1] - 1)
        compute_val_error_path(k, 𝐗_train, 𝐗_val, 𝐚_train, 𝐚_val, 𝐛_train, 𝐛_val, z_min, z_max)
    end).(1:𝑘)
    𝐳s = [paths[i][1] for i = 1:𝑘]
    𝐄s = [paths[i][2] for i = 1:𝑘]
    # merge
    𝐳_k = [z_min]
    𝐄_k = (eltype(𝐄s[begin]))[]
    pointers = ones(Int, 𝑘)
    z_left, z_right = z_min, z_min
    while z_right < z_max
        i_next = argmin((i -> 𝐳s[i][pointers[i] + 1]).(1:𝑘))
        z_right = 𝐳s[i_next][pointers[i_next] + 1]

        push!(𝐳_k, z_right)
        push!(𝐄_k, mean(𝐄s[i][pointers[i]] for i = 1:𝑘))

        z_left = z_right
        pointers[i_next] += 1
    end
    @assert z_right == z_max
    return 𝐳_k, 𝐄_k
end

function compute_Z_CV(k_obs, K, 𝐳, 𝐄)
    z_min, z_max = 𝐳[k_obs][begin], 𝐳[k_obs][end]
    # 各kでの𝐳_kでマージソートしていきながら各区間でE_kobsとE_k′を比較してE_kobs = min E_kとなる区間を求める
    Z_CV = ClosedInterval{eltype(𝐳[k_obs])}[] # init Z_CV
    pointers = Dict(k => 1 for k ∈ K)
    z_left, z_right = z_min, z_min
    while z_right < z_max
        # this algorithm is just like merge sort for 𝐳_k (k ∈ K)
        k_next = K[argmin((k -> 𝐳[k][pointers[k] + 1]).(K))]
        z_right = 𝐳[k_next][pointers[k_next] + 1]
        z_left == z_right && (pointers[k_next] += 1; continue)
        I = z_left..z_right # current interval

        Z_I = ClosedInterval{eltype(𝐳[k_obs])}[] # truncated interval on current I
        for k′ ∈ K
            k_obs == k′ && continue
            Z_I = I ∩ (𝐄[k_obs][pointers[k_obs]] ≤ 𝐄[k′][pointers[k′]])
            isempty(Z_I) && break
        end # Z_I = {z ∈ I | E_{k_obs}(z) = min_k E_k(z)}
        Z_CV = vcat(Z_CV, Z_I)

        z_left = z_right
        pointers[k_next] += 1
    end
    @assert z_right == z_max
    return Z_CV
end

function parametric_SFS_SI_TVS(𝐗, 𝐲_obs, 𝚺, K; ratio=0.5, selection_event=Active())
    𝐴_obs, 𝐬_obs, k_obs = sfs_TVS(𝐗, 𝐲_obs, K; ratio=ratio)
    𝐗_obs = 𝐗[:, 𝐴_obs]
    𝜼 = 𝐗_obs / (𝐗_obs' * 𝐗_obs)
    p_values = Vector{Float64}(undef, length(𝐴_obs)) # selective p-value list
    for j = 1:length(𝐴_obs)
        # parametric representation for 𝐲(z)
        z_obs, σ²_z, 𝐚, 𝐛 = SIforSFS.parametric_representation(𝜼[:, j], 𝐲_obs, 𝚺)

        # compute Z_CV = {z | k(z) = k_obs}
        𝐗_train, 𝐗_val = train_val_split(𝐗; ratio=ratio)
        𝐚_train, 𝐚_val = train_val_split(𝐚; ratio=ratio)
        𝐛_train, 𝐛_val = train_val_split(𝐛; ratio=ratio)
        paths = (k -> compute_val_error_path(k, 𝐗_train, 𝐗_val, 𝐚_train, 𝐚_val, 𝐛_train, 𝐛_val, -abs(z_obs) - 3 * √(σ²_z), abs(z_obs) + 3 * √(σ²_z))).(K)
        𝐳 = Dict(k => paths[i][1] for (i, k) ∈ enumerate(K))
        𝐄 = Dict(k => paths[i][2] for (i, k) ∈ enumerate(K))
        Z_CV = compute_Z_CV(k_obs, K, 𝐳, 𝐄)

        # compute Z_alg = {z | 𝐴(z) = 𝐴_obs}
        Z_alg = truncated_interval(selection_event, 𝐴_obs, 𝐬_obs, 𝐗, 𝐚, 𝐛, k_obs, -abs(z_obs) - 3 * √(σ²_z), abs(z_obs) + 3 * √(σ²_z))

        Z = Z_alg ∩ Z_CV
        𝜋_j = cdf(TruncatedDistribution(Normal(0, √(σ²_z)), Z), z_obs)
        p_values[j] = 2 * min(𝜋_j, 1 - 𝜋_j)
    end
    return (𝐴_obs, p_values)
end

function parametric_SFS_SI_CV(𝐗, 𝐲_obs, 𝚺, K; 𝑘=length(𝐲_obs), selection_event=Active())
    𝐴_obs, 𝐬_obs, k_obs = sfs_CV(𝐗, 𝐲_obs, K; 𝑘=𝑘)
    𝐗_obs = 𝐗[:, 𝐴_obs]
    𝜼 = 𝐗_obs / (𝐗_obs' * 𝐗_obs)
    p_values = Vector{Float64}(undef, length(𝐴_obs)) # selective p-value list
    for j = 1:length(𝐴_obs)
        # parametric representation for 𝐲(z)
        z_obs, σ²_z, 𝐚, 𝐛 = SIforSFS.parametric_representation(𝜼[:, j], 𝐲_obs, 𝚺)

        # compute Z_CV = {z | k(z) = k_obs}
        paths = (k -> compute_cv_path(k, 𝐗, 𝐚, 𝐛, 𝑘, -abs(z_obs) - 3 * √(σ²_z), abs(z_obs) + 3 * √(σ²_z))).(K)
        𝐳 = Dict(k => paths[i][1] for (i, k) ∈ enumerate(K))
        𝐄 = Dict(k => paths[i][2] for (i, k) ∈ enumerate(K))
        Z_CV = compute_Z_CV(k_obs, K, 𝐳, 𝐄)

        # compute Z_alg = {z | 𝐴(z) = 𝐴_obs}
        Z_alg = truncated_interval(selection_event, 𝐴_obs, 𝐬_obs, 𝐗, 𝐚, 𝐛, k_obs, -abs(z_obs) - 3 * √(σ²_z), abs(z_obs) + 3 * √(σ²_z))

        Z = Z_alg ∩ Z_CV
        𝜋_j = cdf(TruncatedDistribution(Normal(0, √(σ²_z)), Z), z_obs)
        p_values[j] = 2 * min(𝜋_j, 1 - 𝜋_j)
    end
    return (𝐴_obs, p_values)
end

function parametric_SFS_CI_TVS(𝐗, 𝐲_obs, 𝚺, K; ratio=0.5, selection_event=Active(), α=0.05)
    𝐴_obs, 𝐬_obs, k_obs = sfs_TVS(𝐗, 𝐲_obs, K; ratio=ratio)
    𝐗_obs = 𝐗[:, 𝐴_obs]
    𝜼 = 𝐗_obs / (𝐗_obs' * 𝐗_obs)
    CIs = Vector{ClosedInterval{Float64}}(undef, length(𝐴_obs)) # selective CI list
    for j = 1:length(𝐴_obs)
        # parametric representation for 𝐲(z)
        z_obs, σ²_z, 𝐚, 𝐛 = SIforSFS.parametric_representation(𝜼[:, j], 𝐲_obs, 𝚺)

        # compute Z_CV = {z | k(z) = k_obs}
        𝐗_train, 𝐗_val = train_val_split(𝐗; ratio=ratio)
        𝐚_train, 𝐚_val = train_val_split(𝐚; ratio=ratio)
        𝐛_train, 𝐛_val = train_val_split(𝐛; ratio=ratio)
        paths = (k -> compute_val_error_path(k, 𝐗_train, 𝐗_val, 𝐚_train, 𝐚_val, 𝐛_train, 𝐛_val, -20 * √(σ²_z), 20 * √(σ²_z))).(K)
        𝐳 = Dict(k => paths[i][1] for (i, k) ∈ enumerate(K))
        𝐄 = Dict(k => paths[i][2] for (i, k) ∈ enumerate(K))
        Z_CV = compute_Z_CV(k_obs, K, 𝐳, 𝐄)

        # compute Z_alg = {z | 𝐴(z) = 𝐴_obs}
        Z_alg = truncated_interval(selection_event, 𝐴_obs, 𝐬_obs, 𝐗, 𝐚, 𝐛, k_obs, -20 * √(σ²_z), 20 * √(σ²_z))

        Z = Z_alg ∩ Z_CV

        CIs[j] = confidence_interval(z_obs, √(σ²_z), Z; α=α)
    end
    return (𝐴_obs, CIs)
end

function parametric_SFS_CI_CV(𝐗, 𝐲_obs, 𝚺, K; 𝑘=length(𝐲_obs), selection_event=Active(), α=0.05)
    𝐴_obs, 𝐬_obs, k_obs = sfs_CV(𝐗, 𝐲_obs, K; 𝑘=𝑘)
    𝐗_obs = 𝐗[:, 𝐴_obs]
    𝜼 = 𝐗_obs / (𝐗_obs' * 𝐗_obs)
    CIs = Vector{ClosedInterval{Float64}}(undef, length(𝐴_obs)) # selective CI list
    for j = 1:length(𝐴_obs)
        # parametric representation for 𝐲(z)
        z_obs, σ²_z, 𝐚, 𝐛 = SIforSFS.parametric_representation(𝜼[:, j], 𝐲_obs, 𝚺)

        # compute Z_CV = {z | k(z) = k_obs}
        paths = (k -> compute_cv_path(k, 𝐗, 𝐚, 𝐛, 𝑘, -20 * √(σ²_z), 20 * √(σ²_z))).(K)
        𝐳 = Dict(k => paths[i][1] for (i, k) ∈ enumerate(K))
        𝐄 = Dict(k => paths[i][2] for (i, k) ∈ enumerate(K))
        Z_CV = compute_Z_CV(k_obs, K, 𝐳, 𝐄)

        # compute Z_alg = {z | 𝐴(z) = 𝐴_obs}
        Z_alg = truncated_interval(selection_event, 𝐴_obs, 𝐬_obs, 𝐗, 𝐚, 𝐛, k_obs, -20 * √(σ²_z), 20 * √(σ²_z))

        Z = Z_alg ∩ Z_CV
        CIs[j] = confidence_interval(z_obs, √(σ²_z), Z; α=α)
    end
    return (𝐴_obs, CIs)
end

end # module ParametricCVforSFS

# # test
# N = 100
# β₀ = 0.0
# β = rand(10)
# K = (1, 3, 5, 7, 9)
# ratio = 0.9
# 𝑘 = 5
# using LinearAlgebra
# 𝚺 = I

# 𝐗, 𝐲_obs = make_dataset(β₀, β, N)
# 𝐗, 𝐲_obs = centering(𝐗, 𝐲_obs)
# 𝐴_obs, _, k_obs = sfs_TVS(𝐗, 𝐲_obs, K; ratio=ratio)
# 𝐴_obs, _, k_obs = sfs_CV(𝐗, 𝐲_obs, K; 𝑘=𝑘)

# 𝐗_obs = 𝐗[:, 𝐴_obs]
# 𝜼 = 𝐗_obs / (𝐗_obs' * 𝐗_obs)
# 𝜼ⱼ = 𝜼[:, 1]
# z_obs = 𝜼ⱼ' * 𝐲_obs
# σ²_z = 𝜼ⱼ' * 𝚺 * 𝜼ⱼ # variance of z_obs
# 𝐛 = 𝚺 * 𝜼ⱼ / σ²_z
# 𝐚 = 𝐲_obs - 𝐛 * 𝜼ⱼ' * 𝐲_obs
# @assert 𝐚 + 𝐛 * z_obs ≈ 𝐲_obs
# # 𝐗_train, 𝐗_val = train_val_split(𝐗; ratio=ratio)
# # 𝐚_train, 𝐚_val = train_val_split(𝐚; ratio=ratio)
# # 𝐛_train, 𝐛_val = train_val_split(𝐛; ratio=ratio)
# # paths = (k -> compute_val_error_path(k, 𝐗_train, 𝐗_val, 𝐚_train, 𝐚_val, 𝐛_train, 𝐛_val, -20 * √(σ²_z), 20 * √(σ²_z))).(1:K) # TVS version
# paths = (k -> compute_cv_path(k, 𝐗, 𝐚, 𝐛, 𝑘, -20 * √(σ²_z), 20 * √(σ²_z))).(K) # CV version
# 𝐳 = Dict(k => paths[i][1] for (i, k) ∈ enumerate(K))
# 𝐄 = Dict(k => paths[i][2] for (i, k) ∈ enumerate(K))
# Z_CV = compute_Z_CV(k_obs, K, 𝐳, 𝐄)
# # create Z_new
# Z_new = eltype(Z_CV)[Z_CV[begin]]
# for I in (@view Z_CV[2:end])
#     if abs(Z_new[end].right - I.left) < 1e-2 * √(σ²_z)
#         Z_new[end] = (Z_new[end].left)..(I.right)
#     else
#         push!(Z_new, I)
#     end
# end
# Z_new
# function val_error(𝐳_k, 𝐄_k)
#     T = length(𝐄_k) # change point size
#     input = Vector{Float64}[]
#     output = Vector{Float64}[]
#     for t in 1:T
#         z_grid = range(𝐳_k[t], 𝐳_k[t + 1], length=10000)
#         input = vcat(input, z_grid)
#         output = vcat(output, 𝐄_k[t].(z_grid))
#     end
#     return input, output
# end
# using Plots
# plot()
# minE, maxE = Ref(Inf), Ref(-Inf)
# for k ∈ K
#     input, output = val_error(𝐳[k], 𝐄[k])
#     minE[] = min(minE[], minimum(output))
#     maxE[] = max(maxE[], maximum(output))
#     plot!(input, output; label="k=$(k)")
# end
# for interval in Z_new # (Z_CV)
#     z_cv = range(interval.left, interval.right; length=1000)
#     plot!(z_cv, zero.(z_cv) .+ (minE[] - 0.01); lw=3, c=:red, label=nothing)
# end
# function point!(x, y; kwargs...) # plot 1 data point
#     scatter!([x], [y]; kwargs...)
# end
# function stem!(x, y_min, y_max; kwargs...) # plot x = x line
#     plot!([x, x], [y_min, y_max]; kwargs...)
# end
# point!(z_obs, minE[] - 0.01; c=:red, label=nothing)
# for interval in Z_new
#     stem!(interval.left, minE[] - 0.01, maxE[] + 0.01; label=nothing, c=:gray, α=0.5, ls=:dashdot)
#     stem!(interval.right, minE[] - 0.01, maxE[] + 0.01; label=nothing, c=:gray, α=0.5, ls=:dashdot)
# end
# using LaTeXStrings
# p = plot!(; xlabel=L"z", ylabel=L"E_k(z)", title=L"k_{\mathrm{obs}}=%$(k_obs)")
# savefig(p, "cv_demo.pdf")


text:
Position(env::GoEnv; args...) = GoPosition(env; args...)
Go(n) = GoEnv(n)

Position(env::GomokuEnv; args...) = GomokuPosition(env; args...)


text:
module TestCSR

include(Pkg.dir("SparseMatrices", "src", "SparseMatrices.jl"))

using FactCheck

I = [1, 2, 1]
J = [1, 1, 3]
V = [1, 1, 2]

Ar = SparseMatrices.SparseMatrixCSR(I, J, V, 2, 3)
A_full = [1 0 2; 1 0 0]

facts("Test SparseMatrixCSR") do

    context("test A_mul_B!") do

        context("CSR*StridedVector") do
            x = [1, 2, 3]
            @fact Ar*x --> [7, 1]
            @fact Base.A_mul_B!(1.0, Ar, x, 0.0, [0, 0]) --> [7, 1]
            @fact Base.A_mul_B!(10.0, Ar, x, 0.0, rand(Int, 2)) --> [70, 10]
            @fact Base.A_mul_B!(10.0, Ar, x, 1.0, [1, 1]) --> [71, 11]
            @fact Base.A_mul_B!(rand(Int, 2), Ar, x) --> [7, 1]

            # check error throwing - Dimension mismatch for wrong size x or out
            @fact_throws DimensionMismatch Ar*[1, 2]
            @fact_throws DimensionMismatch Base.A_mul_B!(1.0, Ar, x, 0.0, [1, 1, 1])
        end

        context("CSR * StridedMatrix") do
            x = [1 2; 3 4; 5 6]
            @fact Ar*x --> [11 14; 1 2]
            @fact Base.A_mul_B!(1.0, Ar, x, 0.0, [0 0; 0 0]) -->  [11 14; 1 2]
            @fact Base.A_mul_B!(10.0, Ar, x, 0.0, rand(Int, 2, 2)) --> [110 140; 10 20]
            @fact Base.A_mul_B!(10.0, Ar, x, 1.0, [1 1; 1 1]) --> [111 141; 11 21]
            @fact Base.A_mul_B!(rand(Int, 2, 2), Ar, x) -->  [11 14; 1 2]

            # check error throwing - Dimension mismatch for wrong size x or out
            @fact_throws DimensionMismatch Ar*rand(Int, 2, 2)
            @fact_throws DimensionMismatch Base.A_mul_B!(1.0, Ar, x, 0.0, rand(Int, 3, 2))
            @fact_throws DimensionMismatch Base.A_mul_B!(1.0, Ar, x, 0.0, rand(Int, 2, 3))
        end

    end
end

end  # module


text:
using Test, PlasmaProperties # This load both the test suite and our MyAwesomePackage

@test 122 + 1 == 123   # All the tests remain to be defined...


text:
module FWFTables

import Tables
import Parsers
import Formatting
import FixedSizeStrings.FixedSizeString

export readbla, Varspec, FWFTable, File, write

"""
    Varspec

Struct to hold all relevant information for one variable.
Fields are

- name: name of the variable
- datatype: Julia type of variable
- startpos: starting position of variabele in record
- length: length in bytes of variable
- decimals: for float-type the number of decimals

Predefined usable types are `FixedSizeString{N}` for string variables,
`Union{Missing,Int64}` and `IntXX` for integer variables, `FloatXX` for real numbers and
`Nothing` for variables that must be skipped. If you want to use another type for input,
then you must define `fw_convert` to convert the byte-string to the given type.

To write to a fixed width file, generator functions for the used types must be available.
Predefined function are defined for `FixedSizeString{N}`, `Union{Missing,Int64}`,
`Float64` and `Nothing`. For other functions the function `makewrite` must be defined.
"""
struct Varspec
    name::String
    datatype::Type
    startpos::Int64
    length::Int64
    decimals::Int64
end

varregex = r"(?i)\s*((?P<name>\w+)\s*(?P<tekst>\".*\"|)\s*:|)\s*(?P<array>ARRAY\[(?P<astart>\d+)\.\.(?P<aeind>\d+)\]\s*OF\s*|)(?P<type>DUMMY|STRING|REAL|INTEGER)\s*\[\s*(?P<length>\d+)(\s*,\s*(?P<decimals>\d+)|)\s*\].*"

convdatatype = Dict(
    "string" => FixedSizeString,
    "integer" => Union{Missing,Int64},
    "dummy" => Nothing,
    "real" => Float64,
)

Base.tryparse(::Type{Int}, ::Nothing) = 0

"""
    readbla(filename)

Read a Blaise specification file for a fixed width ascii file. Returns a vector
of 'Varspec' objects.

A simple example of a Blaise specification file looks like
```
datamodel testbla
 Fields
    var1        : integer[9]
    var2        : string[8]
    dummy[2]
    var3        : Real[8, 2]
    vars        : Array[2010..2050] of STRING[4]
endmodel
```

"""
function readbla(filename)
    spec = Vector{Varspec}()
    regels = open(filename) do io
        readlines(io)
    end
    parsevars = false
    startpos = 1
    for regel in regels
        if !parsevars
            if occursin(r"(?i)fields", regel)
                parsevars = true
            end
            continue
        end
        var = match(varregex, regel)
        if !isnothing(var)
            if !isnothing(var[:name])
                name = lowercase(var[:name])
            else
                name = "dummy"
            end
            len = parse(Int64, var[:length])
            decimals = Base.tryparse(Int64, var[:decimals])
            datatype = convdatatype[lowercase(var[:type])]
            if var[:array] == ""
                push!(spec, Varspec(name, datatype, startpos, len, decimals))
                startpos += len
            else
                # vector of vars
                for varnr = parse(Int64, var[:astart]):parse(Int64, var[a:eind])
                    tmpname = name * "_" * string(varnr)
                    push!(spec, Varspec(tmpname, datatype, startpos, len, decimals))
                    startpos += len
                end
            end
        end
    end
    return spec
end

struct CharVector{N,L} <: AbstractVector{FixedSizeString{N}}
    buffer::Vector{UInt8}
    offset::Int64
    recordlength::Int64
end

Base.IndexStyle(cv::CharVector{N,L}) where {N,L} = IndexLinear()
Base.size(cv::CharVector{N,L}) where {N,L} = (L, )

function Base.getindex(cv::CharVector{N,L}, i::Integer) where {N,L}
    @boundscheck 1 <= i <= L || throw("index $i out of bounds")
    startpos = (i - 1) * cv.recordlength + cv.offset
    endpos = (i - 1) * cv.recordlength + cv.offset + N - 1
    @inbounds s = FixedSizeString{N}(view(cv.buffer, startpos:endpos))
    return s
end

#Base.view(cv::CharVector{N, L}, i::Integer) where {N, L} = cv[i]
#Base.view(cv::CharVector{N, L}, inds::AbstractUnitRange) where {N, L} = view(cv, collect(inds))
#Base.getindex(cv::CharVector{N, L}, inds::AbstractUnitRange) where {N, L} = view(cv, collect(inds))
#Base.getindex(cv::CharVector{N, L}, inds) where {N, L} = view(cv, inds)
#Base.Broadcast.dotview(cv::CharVector{N, L}, inds::AbstractUnitRange) where {N, L} = Base.Broadcast.dotview(cv, collect(inds))

function Base.setindex!(cv::CharVector{N,L}, v, i::Integer) where {N,L}
    startpos = (i - 1) * cv.recordlength + cv.offset
    endpos = (i - 1) * cv.recordlength + cv.offset + N - 1
    cv.buffer[startpos:endpos] = UInt8.(collect(v))
end

function Base.setindex!(cv::CharVector{N,L}, v, inds::AbstractUnitRange) where {N,L}
    for i in inds
        startpos = (i - 1) * cv.recordlength + cv.offset
        endpos = (i - 1) * cv.recordlength + cv.offset + N - 1
        cv.buffer[startpos:endpos] = UInt8.(collect(v))
    end
end

function Base.similar(cv::CharVector{N,L}) where {N,L}
    return Vector{FixedSizeString{N}}(undef, L)
end

function Base.copy(cv::CharVector{N,L}) where {N,L}
    res = Vector{FixedSizeString{N}}(undef, L)
    recordlength = cv.recordlength
    pointer_from = pointer(cv.buffer, 1) + cv.offset - 1
    pointer_to = convert(Ptr{UInt8}, pointer(res, 1))
    for i in 1:L
        unsafe_copyto!(pointer_to, pointer_from, N)
        pointer_from = pointer_from + recordlength
        pointer_to = pointer_to + N
    end
    return res
end

Tables.allocatecolumn(::Type{FixedSizeString{N}}, L) where {N} =
    Vector{FixedSizeString{N}}(undef, L)

function Base.copyto!(
    dest::Vector{FixedSizeString{N}},
    d_o::Integer,
    src::CharVector{N,M},
) where {N,M}
    recordlength = src.recordlength
    pointer_from = pointer(src.buffer, 1) + src.start - 1
    pointer_to = convert(Ptr{UInt8}, pointer(dest, d_o))
    for i = 1:M
        unsafe_copyto!(pointer_to, pointer_from, N)
        pointer_from = pointer_from + recordlength
        pointer_to = pointer_to + N
    end
end

struct FWFTable <: Tables.AbstractColumns
    specs::Vector{Varspec}
    columns::Dict{Symbol,AbstractVector}
end

Tables.istable(::FWFTable) = true
Tables.columnaccess(::FWFTable) = true

specs(t::FWFTable) = getfield(t, :specs)
cols(t::FWFTable) = getfield(t, :columns)

Tables.getcolumn(t::FWFTable, i::Int) = cols(t)[Symbol(specs(t)[i].name)]
Tables.getcolumn(t::FWFTable, nm::Symbol) = cols(t)[nm]
Tables.columnnames(t::FWFTable) = [Symbol(spec.name) for spec in specs(t)]
names(t::FWFTable) = [Symbol(var.name) for var in specs(t)]
types(t::FWFTable) = [var.datatype for var in specs(t)]
Tables.schema(t::FWFTable) = Tables.schema(names(t), types(t))
Base.length(t::FWFTable) = length(cols(t)[Symbol(specs(t)[1].name)])

function createcolumn(::Type{FixedSizeString}, buffer, nrow, startpos, length, recordlength)
    CharVector{length, nrow}(buffer, startpos, recordlength)
end

fw_convert(T::Type{<:Integer}, value) = Parsers.parse(T, value)
fw_convert(::Type{Union{Missing, Int64}}, value) = something(Parsers.tryparse(Int64, value), missing)
fw_convert(::Type{Float64}, value) = something(Parsers.tryparse(Float64, value), NaN64)
fw_convert(::Type{Float32}, value) = something(Parsers.tryparse(Float32, value), NaN32)
fw_convert(::Type{Float16}, value) = something(Parsers.tryparse(Float16, value), NaN32)

function createcolumn(T, buffer, nrow, startpos, length, recordlength)
    column = Vector{T}(undef, nrow)
    pos = startpos
    endpos = startpos + length - 1
    for i = 1:nrow
        setindex!(column, fw_convert(T, view(buffer, pos:endpos)), i)
        pos = pos + recordlength
        endpos = endpos + recordlength
    end
    return column
end

"""
    File(filename::String, blafilename::String)

Read a fixed width file, using the specs from blafile. Returns an object
implementing the Tables interface.

If the file starts with the 3 byte BOM for utf-8, it will be skipped.
The line-ending is determined based on the first line. If the file uses a mixed line
ending, the data will not be read correctly.

# Examples
```julia-repl
julia> using DataFrames, FWFTables
julia> df = DataFrame(FWFTables.File("data.asc", "spec.bla"))
```
"""
function File(filename::String, blafilename::String)
    specs = readbla(blafilename)
    File(filename, specs)
end

"""
    File(filename::String, specs::Vector{Varspec})

Read a fixed width file.
"""
function File(filename::String, specs::Vector{Varspec})
    open(filename) do io
        File(io, specs)
    end
end

"""
    File(io:IO, specs::Vector{Varspec})

Read a fixed width file.
"""
function File(io::IO, specs::Vector{Varspec})
    recordlength = maximum(spec.startpos + spec.length - 1 for spec in specs)
    specselectie = [x for x in specs if x.datatype !== Nothing]
    pos = position(io)
    bom = read(io, 3) == [0xef, 0xbb, 0xbf]
    if !bom
        seek(io, pos)
    end
    buffer = read(io)
    rawlength = recordlength
    while (length(buffer) > recordlength) && (buffer[recordlength+1] in [0x0a, 0x0d])
        recordlength = recordlength + 1
    end
    crlflength = recordlength - rawlength
    # geen regeleinde bij laatste regel opvangen
    if crlflength > 0 && !(buffer[end] in [0x0a, 0x0d])
        nrow = (length(buffer) + crlflength) ÷ recordlength
    else
        nrow = length(buffer) ÷ recordlength
    end
    columns = Dict{Symbol,AbstractVector}()
    Threads.@threads for spec in specselectie
        column = createcolumn(spec.datatype, buffer, nrow, spec.startpos, spec.length, recordlength)
        columns[Symbol(spec.name)] = column
    end
    return FWFTable(specselectie, columns)
end

makewrite(::Type{FixedSizeString}, spec) = (io, val) -> Base.write(io, Ref(val))
function makewrite(::Type{Union{Missing, Int64}}, spec)
    fs = "0>" * string(spec.length) * "d"
    fe = Formatting.FormatSpec(fs)
    spaces = repeat(" ", spec.length)
    return (io, val) -> ismissing(val) ? spaces : Formatting.printfmt(io, fe, val)
end
function makewrite(::Type{Float64}, spec)
    fs = "0>" * string(spec.length) * "." * string(spec.decimals) * "f"
    fe = Formatting.FormatSpec(fs)
    return (io, val) -> Formatting.printfmt(io, fe, val)
end
function makewrite(::Type{Nothing}, spec)
    spaces = repeat(" ", spec.length)
    return io -> Base.write(io, spaces)
end

"""
    write(filename::String, blafilename::String, table)

Write the table to a fixed width ascii file. The blafile contains the specifications of
the file in the Blaise format.

Using the datatype, length and if relevant the decimals from the specification a conversion
function is generated to write a single value using the required width.
For String, Integer, Real and Dummy there are conversion-functions generated by
*makewrite*.  For other datatypes, the required generator must be created. The signature
is (T::Type, spec::Varspec), and the return value must be a function that accepts
(io::IO, value::T) and writes the ascii-representation of the value to io.
"""
function write(filename::String, blafilename::String, table)
    specs::Vector{Varspec} = readbla(blafilename)
    write(filename, specs, table)
end

"""
    write(filename::String, specs::Vector{Varspec}, table)

Write the table to a fixed width ascii file
"""
function write(filename::String, specs::Vector{Varspec}, table)
    open(filename, "w") do io
        write(io, specs, table)
    end
end

"""
    write(io::IO, specs::Vector{Varspec}, table)

Write the table to a fixed width ascii file
"""
function write(io::IO, specs::Vector{Varspec}, table)
    writefcies = [Symbol(spec.name) => makewrite(spec.datatype, spec) for spec in specs]
    for row in Tables.rows(table)
        for (name, writefcie) in writefcies
            if name == :dummy
                writefcie(io)
            else
                writefcie(io, row[name])
            end
        end
        println(io, "")
    end
end

end # module

# vim: ts=4:sw=4:textwidth=89:fileencoding=utf-8

max_stars_repo_name: HenricoWitvliet/FWFTables


text:
#=
Combining MPS-SwapTest method and Quantum Gradient Optimization,
train a MPS circuit with adjustable parameters to recreate target wave function.
=#
export MSQRtrain!, GDescent

"""
    MSQRtrain!(regTar::ArrayReg, MSCircuit::ChainBlock, nTrain::Union{Int64, :auto}; nMeasure::Int64=1, Gmethod::String="Qdiff", GDmethod=("default",0.01), show::Bool=false, useCuYao::Bool=CUDA_ON, ConvTh::Tuple{Float64, Float64}=(5e-4, 1e-3))
    ->
    overlaps::Array{Float64,1}
MSQR training function. This function will change the parameters of differentiable gates in MSCircuit. When set `nTrain = :auto`, trigger the automaic training ieration.
"""
function MSQRtrain!(regTar::ArrayReg, MSCircuit::ChainBlock, nTrain::Union{Int64, Symbol}; nMeasure::Int64=1,
                    Gmethod::Union{String, Tuple{String,Float64}}="Qdiff", GDmethod=("default",0.01), show::Bool=false,
                    useCuYao::Bool=CUDA_ON, ConvTh::Tuple{Float64, Float64}=(5e-4, 1e-3))
    cPar = MSCpar(MSCircuit)
    vBit = cPar.vBit
    rBit = cPar.rBit
    if show
        nBitT = cPar.nBitT
        depth = cPar.depth
        println("\nMSQR Training Parameters:")
        println("nBitT=$(nBitT) vBit=$(vBit) rBit=$(rBit) depth=$(depth)")
        println("nMeasure=$(nMeasure) nTrain=$(nTrain) GDmethod=$(GDmethod)\n")
        println("Initial overlap = $(MStest(regTar, MSCircuit, nMeasure=nMeasure, useCuYao=useCuYao).overlap)")
    end
    reg0 = zero_state((vBit+rBit+1), nbatch=nMeasure)
    regT = repeat(copy(regTar), nMeasure)
    regA = join(reg0, regT)
    useCuYao == true && (regA = regA |> cu)
    if typeof(nTrain) == Int64
        res = train!(nTrain, MSCircuit, Tmethod = MSCircuit->MStest(MSCircuit, regAll=regA),
                     Gmethod=Gmethod, GDmethod=GDmethod, show=show)
    elseif nTrain == :auto
        res = train!(MSCircuit, Tmethod = MSCircuit->MStest(MSCircuit, regAll=regA),
                     Gmethod=Gmethod, GDmethod=GDmethod, show=show, ConvTh=ConvTh)
    end
    res
end

"""
    getDiffs!(circuit::ChainBlock, dGates::Array{QMPS.QDiff, 1}, witnessOp::Yao.Add{<:Any}; trainf::function, showStep::Bool=false, Gmethod::Union{String, Tuple{String,Float64}}="Qdiff", j::Int64=1)
    ->
    tm.overlaps::Float64
Core function of train! to get the overlap in each iteration.
"""
@inline function getDiffs!(circuit::ChainBlock, dGates::Array{QMPS.QDiff, 1}, witnessOp::Yao.Add{<:Any},
                   trainf::Function; GDf::Function=GDescent(("default",0.01)),
                   showStep::Bool=false, Gmethod::Union{String, Tuple{String,Float64}}="Qdiff", j::Int64=1)
    if Gmethod == "Qdiff"
        grads = getQdiff!.(()->trainf(circuit).reg, dGates, Ref(witnessOp))
    elseif Gmethod == "Ndiff"
        grads = getNdiff!.(()->trainf(circuit).reg, dGates, Ref(witnessOp))
    elseif Gmethod[1] == "Ndiff"
        grads = getNdiff!.(()->trainf(circuit).reg, dGates, Ref(witnessOp), δ=Gmethod[2])
    end
    dGatesPar = [parameters(dGates[j])[1] for j=1:length(dGates)]
    dGatesPar = GDf(dGatesPar, grads)
    dispatch!.(dGates, dGatesPar)
    tm = trainf(circuit)
    if  showStep
        println("Training Step $(j), overlap = $(tm.overlap)")
    end
    tm.overlap
end

"""
    train!(nTrain::Int64, circuit::ChainBlock; Tmethod::Function, Gmethod::Union{String, Tuple{String,Float64}}="Qdiff", GDmethod=("default",0.01), show::Bool=false)
    ->
    overlaps::Array{Float64,1}
Function that trains parameters of differentiable gates in the input circuit with designated overlap method
in order to make the circuit generate a similar wave function to the target wave function.
"""
function train!(nTrain::Int64, circuit::ChainBlock; Tmethod::Function,
                Gmethod::Union{String, Tuple{String,Float64}}="Qdiff", GDmethod=("default",0.01), show::Bool=false)
    overlaps = Float64[]
    GD = GDescent(GDmethod)
    tm=Tmethod(circuit)
    witnessOp = put(nqubits(tm.witnessOp), 1=>I2) - tm.witnessOp
    dGates = collect_blocks(QMPS.QDiff, circuit)
    for i=1:nTrain
        showCase = (25*i)%nTrain==0
        overlap = getDiffs!(circuit, dGates, witnessOp, Tmethod, GDf=GD, showStep=(showCase && show), Gmethod=Gmethod, j=i)
        push!(overlaps,overlap)
    end
    overlaps
end
"""
    train!(circuit::ChainBlock; Tmethod::Function, Gmethod::Union{String, Tuple{String,Float64}}="Qdiff", GDmethod=("default",0.01), show::Bool=false, ConvTh::Tuple{Float64, Float64}=(5e-4, 1e-3))
    ->
    overlaps::Array{Float64,1}
Method 2 of `train!`:
\nAutomatically stop training when objective function reaches maximum value.
"""
function train!(circuit::ChainBlock; Tmethod::Function,
                Gmethod::Union{String, Tuple{String,Float64}}="Qdiff", GDmethod=("default",0.01), show::Bool=false, ConvTh::Tuple{Float64, Float64}=(5e-4, 1e-3))
    overlaps = Float64[]
    GD = GDescent(GDmethod)
    tm=Tmethod(circuit)
    witnessOp = put(nqubits(tm.witnessOp), 1=>I2) - tm.witnessOp
    dGates = collect_blocks(QMPS.QDiff, circuit)
    nTrain = 200
    gap = 199
    @label train
    for i=nTrain-gap:nTrain
        showCase = (25*i)%nTrain==0
        overlap = getDiffs!(circuit, dGates, witnessOp, Tmethod, GDf=GD, showStep=(showCase && show), Gmethod=Gmethod, j=i)
        push!(overlaps,overlap)
    end
    gap = 99
    if abs(mean(overlaps[nTrain-199:nTrain-100]) - mean(overlaps[nTrain-99:nTrain])) > ConvTh[1] ||
       abs(mean(overlaps[nTrain-199:nTrain-150]) - mean(overlaps[nTrain-49:nTrain])) > ConvTh[1]*2 ||
       overlaps[nTrain-49:nTrain] |> std |> abs > ConvTh[2]
        nTrain = nTrain+100
        @goto train
    end
    overlaps
end

"""
    GDescent(GDmethod::Union{String, Tuple{String,Float64}, Tuple{String,Float64,Tuple{Float64,Float64}}} = ("default",0.01))
    ->
    reF::Function
Function that returns the aimed Gradient Descent method.
"""
function GDescent(GDmethod::Union{String, Tuple{String,Float64}, Tuple{String,Float64,Tuple{Float64,Float64}}}=("default",0.01))
    if GDmethod == "ADAM"
        resFunc = (dGatesPar, grads) -> Flux.Optimise.update!(ADAM(), dGatesPar, grads)
    elseif GDmethod[1] == "ADAM" && (2 ≤ length(GDmethod) ≤ 3)
        resFunc = (dGatesPar, grads) -> Flux.Optimise.update!(ADAM(GDmethod[2:end]...), dGatesPar, grads)
    elseif GDmethod[1] == "default" && length(GDmethod) == 2
        resFunc = (dGatesPar, grads) -> dGatesPar - grads.*GDmethod[2]
    else
        error("The optimizer of Gradient Descent method is not supported or not specified.")
    end
    resFunc
end