text:
export save, load
using SparseArrays

mutable struct FileReader{T}
    file::AbstractString
    chunktype::Type{T}
    data_offset::Int
    mmap::Bool
end

"""
    save(io::IO, val)

Save a value into the IO buffer. In the case of arrays and sparse
matrices, this will save it in a memory-mappable way.

`load(io::IO, t::Type, domain)` will load the object given its domain
"""
function save(ctx, io::IO, val)
    error("Save method for $(typeof(val)) not defined")
end

###### Save chunks ######

const PARTSPEC = 0x00
const CAT = 0x01

# subparts are saved as Parts

"""
    save(ctx, chunk::Union{Chunk, Thunk}, file_path::AbsractString)

Save a chunk to a file at `file_path`.
"""
function save(ctx, chunk::Union{Chunk, Thunk}, file_path::AbstractString)
    open(file_path, "w") do io
        save(ctx, io, chunk, file_path)
    end
end

"""
special case distmem writing - write to disk on the process with the chunk.
"""
function save(ctx, chunk::Chunk{X,DRef}, file_path::AbstractString) where X
    pid = chunk.handle.where

    remotecall_fetch(pid, file_path, chunk.handle) do path, rref
        open(path, "w") do io
            save(ctx, io, chunk, file_path)
        end
    end
end

function save(ctx, io::IO, chunk::Chunk, file_path)
    meta_io = IOBuffer()

    serialize(meta_io, (chunktype(chunk), domain(chunk)))
    meta = take!(meta_io)

    write(io, PARTSPEC)
    write(io, length(meta))
    write(io, meta)
    data_offset = position(io)

    save(ctx, io, collect(ctx, chunk))

    Chunk(chunktype(chunk), domain(chunk), FileReader(file_path, chunktype(chunk), data_offset, false), false)
end

function save(ctx, io::IO, chunk::DArray, file_path::AbstractString, saved_parts::AbstractArray)

    metadata = (chunktype(chunk), domain(chunk), saved_parts)

    # save yourself
    write(io, CAT)
    serialize(io, metadata)

    DArray(metadata...)
    # write each child
end

function save(ctx, io::IO, chunk::DArray, file_path)
    dir_path = file_path*"_data"
    if !isdir(dir_path)
        mkdir(dir_path)
    end

    # save the chunks
    saved_parts = [save(ctx, c, joinpath(dir_path, lpad(i, 4, "0")))
        for (i, c) in enumerate(chunks(chunk))]

    save(ctx, io, chunk, file_path, saved_parts)
    # write each child
end

function save(ctx, chunk::Chunk{X, FileReader}, file_path::AbstractString) where X
   if abspath(file_path) == abspath(chunk.reader.file)
       chunk
   else
       cp(chunk.reader.file, file_path)
       Chunk(chunktype(chunk), domain(chunk),
          FileReader(file_path, chunktype(chunk),
                     chunk.reader.data_offset, false), false)
   end
end

save(chunk::Union{Chunk, Thunk}, file_path::AbstractString) = save(Context(), chunk, file_path)

###### Load chunks ######

"""
    load(ctx::Context, file_path)

Load an Union{Chunk, Thunk} from a file.
"""
function load(ctx::Context, file_path::AbstractString; mmap=false)

    open(file_path) do f
        part_typ = read(f, UInt8)
        if part_typ == PARTSPEC
            c = load(ctx, Chunk, file_path, mmap, f)
        elseif part_typ == CAT
            c = load(ctx, DArray, file_path, mmap, f)
        else
            error("Could not determine chunk type")
        end
    end
    c
end

"""
    load(ctx::Context, ::Type{Chunk}, fpath, io)

Load a Chunk object from a file, the file path
is required for creating a FileReader object
"""
function load(ctx::Context, ::Type{Chunk}, fname, mmap, io)
    meta_len = read(io, Int)
    io = IOBuffer(read(io, meta_len))

    (T, dmn, sz) = deserialize(io)

    DArray(Chunk(T, dmn, sz,
        FileReader(fname, T, meta_len+1, mmap), false))
end

function load(ctx::Context, ::Type{DArray}, file_path, mmap, io)
    dir_path = file_path*"_data"

    metadata = deserialize(io)
    c = DArray(metadata...)
    for p in chunks(c)
        if isa(p.handle, FileReader)
            p.handle.mmap = mmap
        end
    end
    DArray(c)
end

###### Save and Load for actual data #####

function save(ctx::Context, io::IO, m::Array)
    write(io, reinterpret(UInt8, m, (sizeof(m),)))
    m
end

function save(ctx::Context, io::IO, m::BitArray)
    save(ctx, io, convert(Array{Bool}, m))
end

function collect(ctx::Context, c::Chunk{X,FileReader{T}}) where {X,T<:Array}
    h = c.handle
    io = open(h.file, "r+")
    seek(io, h.data_offset)
    arr = h.mmap ? Mmap.mmap(io, h.chunktype, size(c.domain)) :
        reshape(reinterpret(eltype(T), read(io)), size(c.domain))
    close(io)
    arr
end

function collect(ctx::Context, c::Chunk{X, FileReader{T}}) where {X,T<:BitArray}
    h = c.handle
    io = open(h.file, "r+")
    seek(io, h.data_offset)

    arr = h.mmap ? Mmap.mmap(io, Bool, size(c.domain)) :
        reshape(reinterpret(Bool, read(io)), size(c.domain))
    close(io)
    arr
end

function save(ctx::Context, io::IO, m::SparseMatrixCSC{Tv,Ti}) where {Tv, Ti}
    write(io, m.m)
    write(io, m.n)
    write(io, length(m.nzval))

    typ_io = IOBuffer()
    serialize(typ_io, (Tv, Ti))
    buf = take!(typ_io)
    write(io, sizeof(buf))
    write(io, buf)

    write(io, reinterpret(UInt8, m.colptr, (sizeof(m.colptr),)))
    write(io, reinterpret(UInt8, m.rowval, (sizeof(m.rowval),)))
    write(io, reinterpret(UInt8, m.nzval, (sizeof(m.nzval),)))
    m
end

function collect(ctx::Context, c::Chunk{X, FileReader{T}}) where {X, T<:SparseMatrixCSC}
    h = c.handle
    io = open(h.file, "r+")
    seek(io, h.data_offset)

    m = read(io, Int)
    n = read(io, Int)
    nnz = read(io, Int)

    typ_len = read(io, Int)
    typ_bytes = read(io, typ_len)
    (Tv, Ti) = deserialize(IOBuffer(typ_bytes))

    pos = position(io)
    colptr = Mmap.mmap(io, Vector{Ti}, (n+1,), pos)

    pos += sizeof(Ti)*(n+1)
    rowval = Mmap.mmap(io, Vector{Ti}, (nnz,), pos)

    pos += sizeof(Ti)*nnz
    nnzval = Mmap.mmap(io, Vector{Tv}, (nnz,), pos)
    close(io)

    SparseMatrixCSC(m, n, colptr, rowval, nnzval)
end

function getsub(ctx::Context, c::Chunk{X,FileReader{T}}, d) where {X,T<:AbstractArray}
    Chunk(collect(ctx, c)[d])
end

#### Save computation

struct Save <: Computation
    input
    name::AbstractString
end

function save(p::Computation, name::AbstractString)
    Save(p, name)
end

function stage(ctx::Context, s::Save)
    x = cached_stage(ctx, s.input)
    dir_path = s.name * "_data"
    if !isdir(dir_path)
        mkdir(dir_path)
    end
    function save_part(idx, data)
        p = tochunk(data)
        path = joinpath(dir_path, lpad(idx, 4, "0"))
        saved = save(ctx, p, path)

        # release reference created for the purpose of save
        release_token(p.handle)
        saved
    end

    saved_parts = similar(chunks(x), Thunk)
    for i=1:length(chunks(x))
        saved_parts[i] = Thunk(save_part, i, chunks(x)[i])
    end

    sz = size(chunks(x))
    function save_cat_meta(chunks...)
        f = open(s.name, "w")
        saved_parts = reshape(Union{Chunk, Thunk}[c for c in chunks], sz)
        res = save(ctx, f, x, s.name, saved_parts)
        close(f)
        res
    end

    # The DAG has to block till saving is complete.
    res = Thunk(save_cat_meta, saved_parts...; meta=true)
end


text:
function divider(s, t)
    d = s / t
    return d
end

ext: jl

max_stars_repo_name: wigging/julia-computing

max_issues_repo_name: wigging/julia-computing

max_forks_repo_name: wigging/julia-computing


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

using Markdown
using InteractiveUtils

# ╔═╡ 8e3d151e-0b8b-11eb-1f36-319328a7772c
using Flux

# ╔═╡ 92dfefce-0b8b-11eb-32b4-391aecbbb876
using Distributions

# ╔═╡ 54417ed0-0b8c-11eb-2064-230c0020f268
using IterTools: ncycle

# ╔═╡ 622d190a-0b8c-11eb-1f00-63d9bfdca449
using Random

# ╔═╡ 65e0990a-0b8c-11eb-01ce-233699db160b
using LinearAlgebra

# ╔═╡ a68b270e-0b8c-11eb-050e-b704f7137d1c
using Parameters: @with_kw

# ╔═╡ 20ef48ee-0b8e-11eb-3f94-313efbfc9478
using PlutoUI

# ╔═╡ 055b2642-0b8f-11eb-2777-c1159e9e0e21
using Plots

# ╔═╡ 0f73843e-0b90-11eb-1425-5dce556a16a9
using BenchmarkTools

# ╔═╡ c63890b6-0b8a-11eb-3771-5bde29bc994b
md"""
# Regularization

Weight attenuation is equivalent to regularization of $L_2$ norm, Regularization is a common method to deal with overfitting by adding penalty term to model loss function to make the learned model parameter values smaller.

$l(w, b) + \frac{\lambda}{2n}||w||^2$

among them，$||w||^2 = w_{1}^2 + w_{2}^{2}$

## Higer-dimensional linear regression experiments

Taking high dimensional linear regression as an example, an overfitting problem is introduced.The feature dimension of the data sample is $P$, and the set feature of the data of the training set and the test set is $x_1, x_2, ..., x_p$.

$y = 0.05 + \sum_{i=1}^{p} 0.01x_i + \epsilon$

so we set dimension $p = 200$, smaple number $n=20$:

"""

# ╔═╡ 81ed2efc-0b8b-11eb-1b99-33cb760768fb
md"""
- first import some packages:
"""

# ╔═╡ 6d6b9062-0b8c-11eb-1c2c-d1a65f7ff889
md"""
- do simulation:
"""

# ╔═╡ 949227d2-0b8c-11eb-3ce8-211d2b5010b1
begin
	@with_kw mutable struct Args
		n::Int = 20
		p::Int = 200
		lr::Float64 = 0.003
	end

	args = Args()

	w, b = ones(args.p)*0.01, 0.05
	features = rand(Normal(0,1), (120, args.p))
	labels = features * w .+ b + rand(120)/100
	train_data = Flux.Data.DataLoader((Array(features[1:20, :]'),Array(labels[1:20]')),batchsize=1)
	test_x,test_y = Array(features[21:120, :]'), Array(labels[21:120]')
end

# ╔═╡ 37796f14-0b8d-11eb-350a-cb297a03f349
begin
	train_loss = []
	test_loss = []
	for λ in 0:0.1:4
    	model = Chain(Dense(200,1))
    	ps = Flux.params(model)
    	opt = Descent(args.lr)
    	sqnorm(x) = sum(abs2, x)
    	loss(x,y) = Flux.mse(model(x), y) + λ*sum(sqnorm, Flux.params(model))
    	Flux.train!(loss, ps, train_data, opt);
    	push!(train_loss, Flux.mse(model(Array(features[1:20, :]')), Array(labels[1:20]')))
    	push!(test_loss, Flux.mse(model(test_x), test_y))
	end
	plot(0:0.1:4, [train_loss test_loss], label=["train_loss" "test_loss"], xlabel="λ", ylabel="Loss", title="regularization λ influence the loss")
end

# ╔═╡ 80566366-0b8f-11eb-3ee6-37f6064d46fe
md"""
As we can see, the coefficient of regularization is $\lambda$, if $0<\lambda<1$, we have high loss in test but at train set is great,absolutely, we overfitting train set. when the $\lambda$ become more larger, the test loss start to fall down.

## Dropout
In multilayer perceptron, we describe a multilayer perceptron with a single hidden layer, assuming that number of inputs is 4, and the number of hidden layer unit is 5. So

$h_i = \phi(x_1w_{1i} + x_2w_{2i} + x_3w_{3i}+x_4w_{4i} + b_i)$

among them, $\phi$ is activative function, $x_1,...x_4$ is input value, $w_{1i}, ..., w_{4i}$ is hidden layer parameters.

When use dropout way to this hidden layer, hidden units in this layer will have the chance to be dropped.set the probability is $p$. drop out don't change the expect.

Because in training process, hidden units will be dropped out by $p$, so output calculate can't rely too much on anyone in $h_1,..h_5$.

**For example:**
"""

# ╔═╡ ebc4d248-0b8f-11eb-1123-810b0d1586ab
md"**_Import data:_**"

# ╔═╡ 38c61114-0b90-11eb-3c6f-2d622821ce69
begin
	images = Flux.Data.FashionMNIST.images()
	labels_i = Flux.Data.FashionMNIST.labels()
	k = sort(unique(labels_i))
	onehot_labels = Flux.onehotbatch(labels_i, k)
	train_x = rand(784, 6000)
	for i in 1:6000
		train_x[:,i] = Float64.(reshape(images[i],(784,1)))
	end
	train_data_i = Flux.Data.DataLoader(Array(train_x[:, 1:4000]), Array(onehot_labels[:, 1:4000]), batchsize=10, shuffle=true)
end

# ╔═╡ 8fdaa60c-0b90-11eb-0fca-510bd890f89e
md"**_dropout some layers_**:
every layers has the same probability $p$ of being dropped."

# ╔═╡ d7dc8f6a-0b90-11eb-371b-d9be51947ec4
begin
	m_d = Chain(
    Dense(784, 256, relu),
    Dropout(0.2),
    Dense(256, 256, relu),
    Dropout(0.5),
    Dense(256, 10, relu))
	loss1(x, y) = Flux.logitcrossentropy(m_d(x), y)
	opt_d = Descent(0.05)
	ps_d = Flux.params(m_d)
	@btime Flux.train!(loss1, ps_d, ncycle(train_data_i, 10), opt_d)
end

# ╔═╡ 850218e0-0b91-11eb-0fd9-f13ca40fcb84
md"define accuracy function:"

# ╔═╡ 9f6b4ccc-0b91-11eb-09f4-eb564e81dfc1
accuracy(x, y, model) = Flux.mean(Flux.onecold(model(x)) .== Flux.onecold(y))

# ╔═╡ b9dc1886-0b91-11eb-1bc3-9bb26f25af96
md"show the train accuracy:"

# ╔═╡ aa557fa6-0b91-11eb-0743-233eff7ac496
accuracy(Array(train_x[:, 1:4000]), Array(onehot_labels[:, 1:4000]), m_d)

# ╔═╡ d3021716-0b91-11eb-012c-2f69cb71d076
md"show the test accuracy:"

# ╔═╡ e6b006d0-0b91-11eb-2d8d-29558ea1d719
accuracy(Array(train_x[:, 4000:6000]), Array(onehot_labels[:, 4000:6000]), m_d)

# ╔═╡ Cell order:
# ╟─c63890b6-0b8a-11eb-3771-5bde29bc994b
# ╟─81ed2efc-0b8b-11eb-1b99-33cb760768fb
# ╠═8e3d151e-0b8b-11eb-1f36-319328a7772c
# ╠═92dfefce-0b8b-11eb-32b4-391aecbbb876
# ╠═54417ed0-0b8c-11eb-2064-230c0020f268
# ╠═622d190a-0b8c-11eb-1f00-63d9bfdca449
# ╠═65e0990a-0b8c-11eb-01ce-233699db160b
# ╠═a68b270e-0b8c-11eb-050e-b704f7137d1c
# ╠═20ef48ee-0b8e-11eb-3f94-313efbfc9478
# ╠═055b2642-0b8f-11eb-2777-c1159e9e0e21
# ╠═0f73843e-0b90-11eb-1425-5dce556a16a9
# ╟─6d6b9062-0b8c-11eb-1c2c-d1a65f7ff889
# ╠═949227d2-0b8c-11eb-3ce8-211d2b5010b1
# ╠═37796f14-0b8d-11eb-350a-cb297a03f349
# ╟─80566366-0b8f-11eb-3ee6-37f6064d46fe
# ╟─ebc4d248-0b8f-11eb-1123-810b0d1586ab
# ╠═38c61114-0b90-11eb-3c6f-2d622821ce69
# ╟─8fdaa60c-0b90-11eb-0fca-510bd890f89e
# ╠═d7dc8f6a-0b90-11eb-371b-d9be51947ec4
# ╠═850218e0-0b91-11eb-0fd9-f13ca40fcb84
# ╟─9f6b4ccc-0b91-11eb-09f4-eb564e81dfc1
# ╟─b9dc1886-0b91-11eb-1bc3-9bb26f25af96
# ╠═aa557fa6-0b91-11eb-0743-233eff7ac496
# ╠═d3021716-0b91-11eb-012c-2f69cb71d076
# ╠═e6b006d0-0b91-11eb-2d8d-29558ea1d719

max_stars_repo_name: ZhouZhuofei/juliaDL

max_issues_repo_name: ZhouZhuofei/juliaDL

max_forks_repo_name: ZhouZhuofei/juliaDL


text:
using OrdinaryDiffEq, SteadyStateDiffEq, DiffEqParamEstim, Optim, Test

function f(du,u,p,t)
    α = p[1]
    du[1] = 2 -  α*u[1]
    du[2] = u[1] - 4u[2]
end

p = [2.0]
u0 = zeros(2)
s_prob = SteadyStateProblem(f,u0,p)
s_sol = solve(s_prob,SSRootfind())
s_sol = solve(s_prob,DynamicSS(Tsit5(),abstol=1e-4,reltol=1e-3))

# true data is 1.00, 0.25
data = [1.05, 0.23]
obj = build_loss_objective(s_prob,SSRootfind(),L2Loss([Inf],data),
                           maxiters = Int(1e8),
                           abstol = 1e-10, reltol = 1e-10, verbose=true)
result = Optim.optimize(obj, [2.0], Optim.BFGS())
@test result.minimizer[1] ≈ 2.0 atol=2e-1


text:
## read_bin function with full list of argument

function read_bin(fil::String,kt,kk,prec::DataType);

  if ~isempty(kt);
    error("non-empty kt option not implemented yet");
  end;

  if ~isempty(kk);
    error("non-empty kk option not implemented yet");
  end;

  (n1,n2)=MeshArrays.ioSize;

  if prec==Float64;
    reclen=8;
  else;
    reclen=4;
  end;
  tmp1=stat(fil);
  n3=Int64(tmp1.size/n1/n2/reclen);

  fid = open(fil);
  fld = Array{prec,1}(undef,(n1*n2*n3));
  read!(fid,fld);
  fld = hton.(fld);

  n3>1 ? s=(n1,n2,n3) : s=(n1,n2)
  v0=reshape(fld,s);
  convert2gcmfaces(v0);

end

## read_bin function with reduced list of argument

function read_bin(fil::String,prec::DataType);
  read_bin(fil,[],[],prec);
end

## read_bin function with reduced list of argument

function read_bin(fil::String);
  read_bin(fil,[],[],Float32);
end


text:
type Variable{T}
    " Unique Identifier for the variable "
    id::Int64

    " Name of the Variable "
    name::String

    " Current values for this variable indexed by [qp] "
    value::Array{T,1}

    " Current gradients for this variable indexed by [qp] "
    grad::Array{Vec{2,T},1}

    " Current shape functions for this variable indexed by [qp][i] "
    phi::Array{Array{Float64,1}}

    " Current shape function gradients for this variable indexed by [qp][i] "
    grad_phi::Array{Array{Vec{2,Float64}},1}

    " Determinant of the Jacobian pre-multiplied by the weights "
    JxW::Array{Float64,1}

    " Current DoFs "
    dofs::Array{Int64}

    " Current number of dofs on this element "
    n_dofs::Int64

    " Current number of quadrature points on this element "
    n_qp::Int64

    " The value at a node (for use in NodalBCs) "
    nodal_value::T

    " The dof at a node (for use with NodalBCs) "
    nodal_dof::Int64

    Variable(id::Integer, name::String) = new(id, name,
                                              Array{T}(0),
                                              Array{Vec{2,T}}(0),
                                              Array{Array{Float64}}(0),
                                              Array{Array{Vec{2,Float64}},1}(0),
                                              Array{Float64,1}(0),
                                              Array{Int64}(0),
                                              0,
                                              0,
                                              zero(T),
                                              0)
end

" Recompute all of the data inside of a Variable for a given element "
function reinit!{T}(var::Variable, fe_values::FECellValues, dof_indices::Array{Int64}, dof_values::Array{T})
    var.phi = fe_values.N
    var.grad_phi = fe_values.dNdx
    var.JxW = fe_values.detJdV

    n_qp = length(fe_values.detJdV)

    var.value = [function_value(fe_values, qp, dof_values) for qp in 1:n_qp]
    var.grad = [function_gradient(fe_values, qp, dof_values) for qp in 1:n_qp]

    var.dofs = dof_indices
    var.n_dofs = length(dof_values)
    var.n_qp = n_qp
end

" Recompute all of the data inside of a Variable for a given Node "
function reinit!{T}(var::Variable, dof_index::Int64, dof_value::T)
    var.nodal_dof = dof_index
    var.nodal_value = dof_value
end


text:
using SignatureGB
using Singular
P, (x1, x2, x3, x4, x5, x6)  = Singular.PolynomialRing(Fp(65521), ["x1", "x2", "x3", "x4", "x5", "x6"]);
gb = SignatureGB.f5([x1, x2, x3, x4, x5, x6]);F = [841*x1^2*x2^2+2146*x1^2*x2*x3+1369*x1^2*x3^2+8281*x1^2*x4^2-14742*x1^2*x4*x5-11830*x1^2*x4*x6+6561*x1^2*x5^2+10530*x1^2*x5*x6+4225*x1^2*x6^2+2184*x1*x2^2*x4-1654*x1*x2^2*x5-1560*x1*x2^2*x6+370*x1*x2*x3*x5-2088*x1*x2*x6^2-14196*x1*x3*x4^2+12636*x1*x3*x4*x5+10140*x1*x3*x4*x6-2664*x1*x3*x6^2-910*x1*x4^2*x6+810*x1*x4*x5*x6+650*x1*x4*x6^2+144*x2^4-1872*x2^2*x3*x4-120*x2^2*x4*x6+25*x2^2*x5^2-360*x2*x5*x6^2+6084*x3^2*x4^2+780*x3*x4^2*x6+25*x4^2*x6^2+1296*x6^4-3886*x1*x2*x4+3364*x1*x2*x6-4958*x1*x3*x4+4292*x1*x3*x6-670*x2*x4*x5+580*x2*x5*x6+4824*x4*x6^2-4176*x6^3+4489*x4^2-7772*x4*x6+3364*x6^2, 1682*x1^2*x2+2146*x1^2*x3+4368*x1*x2*x4-3308*x1*x2*x5-3120*x1*x2*x6+370*x1*x3*x5-2088*x1*x6^2+576*x2^3-3744*x2*x3*x4-240*x2*x4*x6+50*x2*x5^2-360*x5*x6^2-3886*x1*x4+3364*x1*x6-670*x4*x5+580*x5*x6, 2146*x1^2*x2+2738*x1^2*x3+370*x1*x2*x5-14196*x1*x4^2+12636*x1*x4*x5+10140*x1*x4*x6-2664*x1*x6^2-1872*x2^2*x4+12168*x3*x4^2+780*x4^2*x6-4958*x1*x4+4292*x1*x6, 16562*x1^2*x4-14742*x1^2*x5-11830*x1^2*x6+2184*x1*x2^2-28392*x1*x3*x4+12636*x1*x3*x5+10140*x1*x3*x6-1820*x1*x4*x6+810*x1*x5*x6+650*x1*x6^2-1872*x2^2*x3-120*x2^2*x6+12168*x3^2*x4+1560*x3*x4*x6+50*x4*x6^2-3886*x1*x2-4958*x1*x3-670*x2*x5+4824*x6^2+8978*x4-7772*x6, -14742*x1^2*x4+13122*x1^2*x5+10530*x1^2*x6-1654*x1*x2^2+370*x1*x2*x3+12636*x1*x3*x4+810*x1*x4*x6+50*x2^2*x5-360*x2*x6^2-670*x2*x4+580*x2*x6, -11830*x1^2*x4+10530*x1^2*x5+8450*x1^2*x6-1560*x1*x2^2-4176*x1*x2*x6+10140*x1*x3*x4-5328*x1*x3*x6-910*x1*x4^2+810*x1*x4*x5+1300*x1*x4*x6-120*x2^2*x4-720*x2*x5*x6+780*x3*x4^2+50*x4^2*x6+5184*x6^3+3364*x1*x2+4292*x1*x3+580*x2*x5+9648*x4*x6-12528*x6^2-7772*x4+6728*x6];
sort!(F, by = p -> leading_monomial(p))
println(@elapsed SignatureGB.f5(F, verbose=true));exit()


text:
using Oceananigans
using Oceananigans.Units
using Oceananigans.ImmersedBoundaries: ImmersedBoundaryGrid, GridFittedBottom
using Oceananigans.TurbulenceClosures: VerticallyImplicit, Vertical, Horizontal
using Printf
using GLMakie
using JLD2

grid = RectilinearGrid(topology = (Periodic, Bounded, Bounded),
                       size = (128, 128, 16),
                       x = (-500kilometers, 500kilometers),
                       y = (-500kilometers, 500kilometers),
                       z = (-1kilometers, 0),
                       halo = (3, 3, 3))

const Lz = grid.Lz

# Uncomment to put a bump in the grid:
#=
const width = 50kilometers
bump(x, y) = - Lz * (1 - 0.5 * exp(-x^2 / 2width^2))
grid = ImmersedBoundaryGrid(grid, GridFittedBottom(bump))
=#

Δx = grid.Lx / grid.Nx
κh = Δx^4 / 12hours

diffusive_closure = ScalarDiffusivity(ν = 1e-2, κ = 1e-2,
                                      isotropy = Vertical(),
                                      time_discretization = VerticallyImplicit())

horizontal_closure = ScalarBiharmonicDiffusivity(ν = κh, κ = κh, isotropy = Horizontal())

model = HydrostaticFreeSurfaceModel(grid = grid,
                                    coriolis = BetaPlane(latitude = -45),
                                    buoyancy = BuoyancyTracer(),
                                    closure = (diffusive_closure, horizontal_closure),
                                    tracers = (:b, :c),
                                    momentum_advection = WENO5(),
                                    tracer_advection = WENO5(),
                                    free_surface = ImplicitFreeSurface())

# Initial condition: a baroclinically unstable situation!
ramp(y, δy) = min(max(0, y/δy + 1/2), 1)

# Parameters
N² = 4e-6 # [s⁻²] buoyancy frequency / stratification
M² = 8e-8 # [s⁻²] horizontal buoyancy gradient

δy = 50kilometers
δz = 100
Lz = grid.Lz

δc = 2δy
δb = δy * M²
ϵb = 1e-2 * δb # noise amplitude

bᵢ(x, y, z) = N² * z + δb * ramp(y, δy) + ϵb * randn()
cᵢ(x, y, z) = exp(-y^2 / 2δc^2) * exp(-(z + Lz/4)^2 / 2δz^2)

set!(model, b=bᵢ, c=cᵢ)

simulation = Simulation(model, Δt=20minutes, stop_time=30days)

wizard = TimeStepWizard(cfl=0.2, max_change=1.1, max_Δt=simulation.Δt)
simulation.callbacks[:wizard] = Callback(wizard, IterationInterval(1))

function print_progress(sim)
    @printf("[%05.2f%%] i: %d, t: %s, wall time: %s, max(u): (%6.3e, %6.3e, %6.3e) m/s, next Δt: %s\n",
            100 * (sim.model.clock.time / sim.stop_time),
            sim.model.clock.iteration,
            prettytime(sim.model.clock.time),
            prettytime(sim.run_wall_time),
            maximum(abs, sim.model.velocities.u),
            maximum(abs, sim.model.velocities.v),
            maximum(abs, sim.model.velocities.w),
            prettytime(sim.Δt))

    return nothing
end

simulation.callbacks[:print_progress] = Callback(print_progress, IterationInterval(20))

b, c = model.tracers
u, v, w = model.velocities
ζ = Field(∂x(v) - ∂y(u))

simulation.output_writers[:fields] = JLD2OutputWriter(model, (; ζ, b, c),
                                                      schedule = TimeInterval(1hour),
                                                      field_slicer = FieldSlicer(k=grid.Nz),
                                                      prefix = "baroclinic_adjustment",
                                                      force = true)

run!(simulation)

file = jldopen("baroclinic_adjustment.jld2")
iterations = parse.(Int, keys(file["timeseries/t"]))
ζ = [file["timeseries/ζ/$i"][:, :, 1] for i in iterations]
b = [file["timeseries/b/$i"][:, :, 1] for i in iterations]
c = [file["timeseries/c/$i"][:, :, 1] for i in iterations]
t = [file["timeseries/t/$i"] for i in iterations]
close(file)

Nt = length(t)

fig = Figure()

axζ = Axis(fig[1, 1])
axb = Axis(fig[1, 2])
axc = Axis(fig[1, 3])

slider = Slider(fig[2, :], range=1:Nt, startvalue=1)
n = slider.value

xζ, yζ, zζ = 1e3 .* nodes((Face, Face, Center), grid)
xc, yc, zc = 1e3 .* nodes((Center, Center, Center), grid)

ζⁿ = @lift ζ[$n]
bⁿ = @lift b[$n]
cⁿ = @lift c[$n]

hmζ = heatmap!(axζ, xζ, yζ, ζⁿ)
hmb = heatmap!(axb, xc, yc, bⁿ)
hmc = heatmap!(axc, xc, yc, cⁿ)

# Colorbar(fig[2, 1], xζ, yζ, hmζ, horizontal=true)
# Colorbar(fig[2, 2], xc, yc, hmb, horizontal=true)
# Colorbar(fig[2, 3], xc, yc, hmc, horizontal=true)

display(fig)


text:
function MPI.RBuffer(senddata::ITensor, recvdata::ITensor)
  senddata = permute(senddata, inds(recvdata); allow_alias=true)
  @assert inds(senddata) == inds(recvdata)
  _senddata = ITensors.data(senddata)
  _recvdata = ITensors.data(recvdata)
  return MPI.RBuffer(_senddata, _recvdata)
end

struct MPISum{T}
  data::T
  comm::MPI.Comm
end

MPISum(data::T, comm=MPI.COMM_WORLD) where {T} = MPISum{T}(data, comm)

nsite(P::MPISum) = nsite(P.data)

Base.length(P::MPISum) = length(P.data)

function _Allreduce(sendbuf::ITensor, op, comm)
  return _Allreduce(NDTensors.storagetype(tensor(sendbuf)), sendbuf, op, comm)
end

function _Allreduce(::Type{<:BlockSparse}, sendbuf::ITensor, op, comm)
  # Similar to:
  #
  # recvbuf = ITensor(eltype(sendbuf), flux(sendbuf), inds(sendbuf))
  #
  # But more efficient.
  nprocs = MPI.Comm_size(comm)
  N = order(sendbuf)
  max_nblocks = maximum(MPI.Allgather(nnzblocks(sendbuf), comm))
  blocks_sendbuf = fill(0, max_nblocks, N)
  nzblocks_matrix = Matrix{Int}(reduce(hcat, collect.(nzblocks(sendbuf)))')
  blocks_sendbuf[1:size(nzblocks_matrix, 1), :] = nzblocks_matrix
  allblocks = MPI.Allgather(blocks_sendbuf, comm)
  allblocks_tensor = reshape(allblocks, max_nblocks, N, nprocs)
  allblocks_vector = [allblocks_tensor[:, :, i] for i in 1:nprocs]
  allblocks_vector_tuple = [
    [Tuple(allblocks_vector[n][i, :]) for i in 1:max_nblocks] for n in 1:nprocs
  ]
  blocks = Block.(filter(≠(ntuple(Returns(0), N)), union(allblocks_vector_tuple...)))
  recvbuf = itensor(BlockSparseTensor(eltype(sendbuf), blocks, inds(sendbuf)))
  sendbuf_filled = copy(recvbuf)
  sendbuf_filled .= sendbuf
  MPI.Allreduce!(sendbuf_filled, recvbuf, +, comm)
  return recvbuf
end

function _Allreduce(::Type{<:Dense}, sendbuf::ITensor, op, comm)
  recvbuf = similar(sendbuf)
  MPI.Allreduce!(sendbuf, recvbuf, op, comm)
  return recvbuf
end

function product(P::MPISum, v::ITensor)
  return _Allreduce(P.data(v), +, P.comm)
end

function Base.eltype(P::MPISum)
  return eltype(P.data)
end

(P::MPISum)(v::ITensor) = product(P, v)

Base.size(P::MPISum) = size(P.data)

function position!(P::MPISum, psi::MPS, pos::Int)
  position!(P.data, psi, pos)
  return P
end

function noiseterm(P::MPISum, phi::ITensor, dir::String)
  return _Allreduce(noiseterm(P.data, phi, dir), +, P.comm)
end


text:
using NetworkInference
using Test
using DelimitedFiles

# Only run tests if the EmpiricalBayes package exists
if NetworkInference.EB_EXISTS

using EmpiricalBayes

println("Testing empirical Bayes glue functions...")

data_folder_path = joinpath(dirname(@__FILE__), "data")

# Test to_index functions
@test to_index("bbb", "aaa") == to_index("aaa", "bbb")
@test to_index("aaa", "bbb") == ("bbb", "aaa")

n1 = Node("aaa", [], 0, [])
n2 = Node("bbb", [], 0, [])
@test to_index(n1, n2) == to_index(n2, n1)
@test to_index(n1, n2) == ("bbb", "aaa")

@test to_index([n1, n2]) == to_index([n2, n1])
@test to_index([n1, n2]) == ("bbb", "aaa")

# Test the make_priors function
prior_path = joinpath(data_folder_path, "test_priors.txt")
reference_priors = Dict( [ (("bbb", "aaa"), 1) ,
                     (("ccc", "aaa"), 0) ,
                     (("ccc", "bbb"), 1)
                   ])
@test make_priors(prior_path) == reference_priors
reference_priors = Dict( [ (("ccc", "aaa"), 0) ,
                     (("ccc", "bbb"), 1) ,
                     (("bbb", "aaa"), 0)
                   ])
@test make_priors(prior_path, 1.0, 2) == reference_priors
reference_priors = Dict( [ (("bbb", "aaa"), 2) ,
                     (("ccc", "aaa"), 0) ,
                     (("ccc", "bbb"), 2)
                   ])
@test make_priors(prior_path, 2.0) == reference_priors
@test make_priors([prior_path, prior_path]) == reference_priors
@test make_priors([prior_path, prior_path], [1.5, 0.5]) == reference_priors
reference_priors = Dict( [ (("ccc", "aaa"), 0) ,
                     (("ccc", "bbb"), 2) ,
                     (("bbb", "aaa"), 0)
                   ])
@test make_priors([prior_path, prior_path], [1.5, 0.5], 2) == reference_priors

# Test the empirical_bayes function
# Test with w0 = 2.2, because that is the default for empirical_bayes within NetworkInference
println("Inferring test empirical Bayes networks...")
mi_benchmark = readdlm(joinpath(data_folder_path, "mi.txt"))
mi_benchmark = mi_benchmark[1:2:end, :] # skip repeated edges

mi_priors_filepath = joinpath(data_folder_path, "mi_priors.txt")
mi_priors = readdlm(mi_priors_filepath)
mi_priors = mi_priors[1:2:end, :] # skip repeated edges
prior_dict = make_priors(mi_priors_filepath)

yeast_test_data = joinpath(data_folder_path, "yeast1_10_data.txt")
nodes = get_nodes(yeast_test_data)
mi_network = InferredNetwork(MINetworkInference(), nodes)

eb_network = empirical_bayes(mi_network, prior_dict, 5, :Gamma, tail = :two, w0 = 2.2)
eb_weights = [e.weight for e in eb_network.edges]

benchmark_stats = convert(Array{Float64}, mi_benchmark[:, 3])
benchmark_priors = convert(Array{Float64}, mi_priors[:, 3])
ref_weights = empirical_bayes(benchmark_stats, benchmark_priors, 5, :Gamma, tail = :two, w0 = 2.2)

@test eb_weights ≈ sort(ref_weights, rev = true) atol = 0.0001

# Test the empirical_bayes function with no priors
println("Inferring test empirical Bayes networks with no priors...")
eb_no_prior_network = empirical_bayes(mi_network, 5, :Gamma, tail = :two, w0 = 2.2)
eb_no_prior_weights = [e.weight for e in eb_no_prior_network.edges]

ref_no_prior_weights = empirical_bayes(benchmark_stats, 5, :Gamma, tail = :two, w0 = 2.2)

@test eb_no_prior_weights ≈ sort(ref_no_prior_weights, rev = true) atol = 0.0001

println("Empirical Bayes glue tests passed")

end


text:
module UnitTests

using Distributed
using Base.Filesystem
@everywhere using Pkg
@everywhere Pkg.activate(joinpath(@__DIR__, "..", "..", ".."))
@everywhere Pkg.instantiate()
@everywhere using XUnit

function run_unittests()
    local tests_exception = nothing
    prev_TESTSET_PRINT_ENABLE = XUnit.TESTSET_PRINT_ENABLE[]
    XUnit.TESTSET_PRINT_ENABLE[] = false
    try
        XUnit.runtests(joinpath(@__DIR__, "unittests.jl"))
    catch e
        # swallow the error
        # @error "XUnit Tests Progress" exception=e
    finally
        XUnit.TESTSET_PRINT_ENABLE[] = prev_TESTSET_PRINT_ENABLE
    end

    @testset "unittests tests" begin
        @test isdefined(Main, :XUNIT_UNITTEST_RESULTS)

        test_results = Main.XUNIT_UNITTEST_RESULTS[]

        @test test_results !== nothing

        @test test_results isa XUnit.AsyncTestSuite

        passes, fails, errors, broken, c_passes, c_fails, c_errors, c_broken =
            XUnit._swallow_all_outputs() do
                XUnit.get_test_counts(test_results)
            end
        # There are no direct tests under the top-level test-suite. That's why
        # `passes`, `fails`, `errors` and `broken` are `0`
        @test passes == 0
        @test fails == 0
        @test errors == 0
        @test broken == 0
        # The following tests check the number of test `passes`, `fails`, `errors` and
        # `broken` in "unittests.jl"
        @test c_passes == 55
        @test c_fails == 11
        @test c_errors == 4
        @test c_broken == 0
    end

    nothing
end

run_unittests()

end


text:
"""
    FixedString{N}(str)

An ascii string allocated in an immutable fixed length buffer of `N` `UInt8`
code points. This type is intended to be binary compatible with C buffers such
as might be declared with `char buf[16]`.

This is particularly useful for interacting with C libraries or serialized data
structures with fixed length strings.

Currently an ascii encoding is assumed.
"""
struct FixedString{N} <: AbstractString
    data::NTuple{N,UInt8}
end

# constructor / convert
function FixedString{N}(s::AbstractString) where N
    if !isascii(s)
        throw(ArgumentError("String $(repr(s)) is not ascii"))
    elseif sizeof(s) > N-1
        throw(ArgumentError("String $(repr(s)) does not fit into $N-1 ascii characters"))
    end
    FixedString{N}(ntuple(i->UInt8(i <= sizeof(s) ? s[i] : '\0'), N))
end

Base.convert(::Type{FixedString{N}}, s::FixedString{N}) where {N} = s
Base.convert(::Type{FixedString{N}}, s::AbstractString) where {N} = FixedString{N}(s)

function Base.convert(::Type{String}, s::FixedString)
    vec_data = UInt8[s.data...]
    String(vec_data[1:ncodeunits(s)])
end

# IO
Base.show(io::IO,  s::FixedString) = print(io, typeof(s), "(", repr(convert(String, s)), ")")

Base.write(io::IO, s::FixedString) = write(io, Ref(s))

function Base.read(io::IO, ::Type{<:FixedString{N}}) where N
    read!(io, Ref{FixedString{N}}())[]
end

# AbstractString interface. It's not quite clear what this is yet...
# https://discourse.julialang.org/t/what-is-the-interface-of-abstractstring/8937

# Plenty of assumptions about ascii chars follow!

function Base.ncodeunits(s::FixedString)
    firstnull = findfirst(isequal(0), s.data)
    firstnull != nothing ? firstnull-1 : sizeof(s.data)
end

Base.codeunit(s::FixedString) = UInt8
@inline function Base.codeunit(s::FixedString, i::Integer)
    @boundscheck checkbounds(s, i)
    s.data[i]
end
Base.isvalid(s::FixedString, i::Integer) = true
#Base.@propagate_inbounds Base.nextind(s::FixedString, i::Int) = i+1
Base.length(s::FixedString) = ncodeunits(s)

function Base.iterate(s::FixedString, i::Integer=0)
    i = nextind(s, i)
    (i <= length(s.data) && s.data[i] != 0x00) || return nothing
    return (Char(s.data[i]), i)
end

max_stars_repo_name: TRIImaging/MagneticResonanceSignals

max_issues_repo_name: TRIImaging/MagneticResonanceSignals

max_forks_repo_name: TRIImaging/MagneticResonanceSignals


text:
########################
# Generic approximation
########################

"""
The `approximation` function returns an operator that can be used to approximate
a function in the function set. This operator maps a grid to a set of coefficients.
"""
approximation(Φ; discrete = true, options...) =
  discrete ?
    discrete_approximation(Φ; options...) :
    continuous_approximation(Φ; options...)

function discrete_approximation(Φ::Dictionary; options...)
    if isbasis(Φ) && hasinterpolationgrid(Φ)
        interpolation(Φ; options...)
    else
        leastsquares(Φ; options...)
    end
end

continuous_approximation(Φ::Dictionary; options...) = inv(gram(Φ; options...))

# Automatically sample a function if an operator is applied to it
function (*)(op::DictionaryOperator, f::Function)
    op * project(src(op), f)
end

project(Φ::GridBasis, f::Function; options...) = sample(Φ, f)

approximate(Φ::Dictionary, f::Function; options...) =
    Expansion(Φ, approximation(Φ; options...) * f)

"""
The 2-argument approximation exists to allow you to transform any
Dictionary coefficients to any other Dictionary coefficients, including
Discrete Grid Spaces.
If a transform exists, the transform is used.
"""
function approximation(A::Dictionary, B::Dictionary, options...)
    if hastransform(A, B)
        return transform(A, B; options...)
    else
        error("Don't know a transformation from ", A, " to ", B)
    end
end

#################
# Interpolation
#################

function interpolation_matrix(Φ::Dictionary, pts, T=coefficienttype(Φ))
    @assert length(Φ) == length(pts)
    evaluation_matrix(Φ, pts, T)
end

interpolation(Φ::Dictionary, args...; options...) =
    interpolation(operatoreltype(Φ), Φ, args...; options...)

interpolation(::Type{T}, Φ::Dictionary; options...) where {T} =
    interpolation(T, Φ, interpolation_grid(Φ); options...)

interpolation(::Type{T}, Φ::Dictionary, grid::AbstractGrid; options...) where {T} =
    interpolation(T, Φ, GridBasis{T}(grid); options...)

function interpolation(::Type{T}, Φ::Dictionary, gb::GridBasis; options...) where {T}
    if hasinterpolationgrid(Φ) && interpolation_grid(Φ) == grid(gb) && hastransform(Φ, gb)
        transform(T, gb, Φ; options...)
    else
        default_interpolation(Φ, gb; options...)
    end
end

default_interpolation(Φ::Dictionary, gb::GridBasis; options...) =
    default_interpolation(operatoreltype(Φ, gb), Φ, gb; options...)

default_interpolation(::Type{T}, Φ::Dictionary, gb::GridBasis; options...) where {T} =
    QR_solver(evaluation(T, Φ, grid(gb)))

interpolate(Φ::Dictionary, pts, f, T=coefficienttype(s)) = Expansion(interpolation(T, Φ, pts) * f)

########################
# Generic least squares
########################

function leastsquares_matrix(Φ::Dictionary, pts, T=coefficienttype(Φ))
    @assert length(Φ) <= length(pts)
    evaluation_matrix(Φ, pts, T)
end

leastsquares(Φ::Dictionary, args...; options...) =
    leastsquares(operatoreltype(Φ), Φ, args...; options...)

function leastsquares(::Type{T}, Φ::Dictionary; samplingfactor = 2, options...) where {T}
    if hasinterpolationgrid(Φ)
        Φ2 = resize(Φ, samplingfactor*length(Φ))
        ls_grid = interpolation_grid(Φ2)
    else
        ls_grid = EquispacedGrid(samplingfactor*length(Φ), support(Φ))
    end
    leastsquares(T, Φ, ls_grid; options...)
end

leastsquares(::Type{T}, Φ::Dictionary, grid::AbstractGrid; options...) where {T} =
    leastsquares(T, Φ, GridBasis{T}(grid); options...)

function leastsquares(T, Φ::Dictionary, gb::GridBasis; options...)
    if hasinterpolationgrid(Φ)
        larger_Φ = resize(Φ, size(gb))
        if interpolation_grid(larger_Φ) == grid(gb) && hastransform(larger_Φ, gb)
            R = restriction(T, larger_Φ, Φ; options...)
            E = transform(T, gb, larger_Φ; options...)
            R * E
        else
            default_leastsquares(Φ, gb, T; options...)
        end
    else
        default_leastsquares(Φ, gb, T; options...)
    end
end

default_leastsquares(Φ::Dictionary, gb::GridBasis, T=operatoreltype(Φ,gb); options...) =
    QR_solver(ArrayOperator(leastsquares_matrix(Φ, grid(gb), T), Φ, gb))


text:
abstract type StructOfConstraints <: MOI.ModelLike end

function _throw_if_cannot_delete(model::StructOfConstraints, vis, fast_in_vis)
    broadcastcall(model) do constrs
        return _throw_if_cannot_delete(constrs, vis, fast_in_vis)
    end
end
function _deleted_constraints(
    callback::Function,
    model::StructOfConstraints,
    vi,
)
    broadcastcall(model) do constrs
        return _deleted_constraints(callback, constrs, vi)
    end
end

function MOI.add_constraint(
    model::StructOfConstraints,
    func::MOI.AbstractFunction,
    set::MOI.AbstractSet,
)
    if MOI.supports_constraint(model, typeof(func), typeof(set))
        return MOI.add_constraint(
            constraints(model, typeof(func), typeof(set)),
            func,
            set,
        )
    else
        throw(MOI.UnsupportedConstraint{typeof(func),typeof(set)}())
    end
end

function constraints(
    model::StructOfConstraints,
    ci::MOI.ConstraintIndex{F,S},
) where {F,S}
    if !MOI.supports_constraint(model, F, S)
        throw(MOI.InvalidIndex(ci))
    end
    return constraints(model, F, S)
end
function MOI.get(
    model::StructOfConstraints,
    attr::Union{MOI.ConstraintFunction,MOI.ConstraintSet},
    ci::MOI.ConstraintIndex,
)
    return MOI.get(constraints(model, ci), attr, ci)
end

function MOI.delete(model::StructOfConstraints, ci::MOI.ConstraintIndex)
    return MOI.delete(constraints(model, ci), ci)
end

function MOI.is_valid(
    model::StructOfConstraints,
    ci::MOI.ConstraintIndex{F,S},
) where {F,S}
    if MOI.supports_constraint(model, F, S)
        return MOI.is_valid(constraints(model, ci), ci)
    else
        return false
    end
end

function MOI.modify(
    model::StructOfConstraints,
    ci::MOI.ConstraintIndex,
    change::MOI.AbstractFunctionModification,
)
    return MOI.modify(constraints(model, ci), ci, change)
end

function MOI.set(
    model::StructOfConstraints,
    attr::Union{MOI.ConstraintFunction,MOI.ConstraintSet},
    ci::MOI.ConstraintIndex,
    func_or_set,
)
    return MOI.set(constraints(model, ci), attr, ci, func_or_set)
end

function MOI.get(
    model::StructOfConstraints,
    loc::MOI.ListOfConstraints,
) where {T}
    return broadcastvcat(model) do v
        return MOI.get(v, loc)
    end
end

function MOI.get(
    model::StructOfConstraints,
    noc::MOI.NumberOfConstraints{F,S},
) where {F,S}
    if MOI.supports_constraint(model, F, S)
        return MOI.get(constraints(model, F, S), noc)
    else
        return 0
    end
end

function MOI.get(
    model::StructOfConstraints,
    loc::MOI.ListOfConstraintIndices{F,S},
) where {F,S}
    if MOI.supports_constraint(model, F, S)
        return MOI.get(constraints(model, F, S), loc)
    else
        return MOI.ConstraintIndex{F,S}[]
    end
end

function MOI.is_empty(model::StructOfConstraints)
    return mapreduce_constraints(MOI.is_empty, &, model, true)
end
function MOI.empty!(model::StructOfConstraints)
    return broadcastcall(MOI.empty!, model)
end

# Can be used to access constraints of a model
"""
broadcastcall(f::Function, model::AbstractModel)

Calls `f(contrs)` for every vector `constrs::Vector{ConstraintIndex{F, S}, F, S}` of the model.

# Examples

To add all constraints of the model to a solver `solver`, one can do
```julia
_addcon(solver, ci, f, s) = MOI.add_constraint(solver, f, s)
function _addcon(solver, constrs::Vector)
    for constr in constrs
        _addcon(solver, constr...)
    end
end
MOIU.broadcastcall(constrs -> _addcon(solver, constrs), model)
```
"""
function broadcastcall end

"""
broadcastvcat(f::Function, model::AbstractModel)

Calls `f(contrs)` for every vector `constrs::Vector{ConstraintIndex{F, S}, F, S}` of the model and concatenate the results with `vcat` (this is used internally for `ListOfConstraints`).

# Examples

To get the list of all functions:
```julia
_getfun(ci, f, s) = f
_getfun(cindices::Tuple) = _getfun(cindices...)
_getfuns(constrs::Vector) = _getfun.(constrs)
MOIU.broadcastvcat(_getfuns, model)
```
"""
function broadcastvcat end

function mapreduce_constraints end

# Macro code

abstract type SymbolFS end
struct SymbolFun <: SymbolFS
    s::Union{Symbol,Expr}
    typed::Bool
end
struct SymbolSet <: SymbolFS
    s::Union{Symbol,Expr}
    typed::Bool
end

# QuoteNode prevents s from being interpolated and keeps it as a symbol
# Expr(:., MOI, s) would be MOI.s
# Expr(:., MOI, $s) would be Expr(:., MOI, EqualTo)
# Expr(:., MOI, :($s)) would be Expr(:., MOI, :EqualTo)
# Expr(:., MOI, :($(QuoteNode(s)))) is Expr(:., MOI, :(:EqualTo)) <- what we want

# (MOI, :Zeros) -> :(MOI.Zeros)
# (:Zeros) -> :(MOI.Zeros)
_set(s::SymbolSet) = esc(s.s)
_fun(s::SymbolFun) = esc(s.s)
function _typedset(s::SymbolSet)
    if s.typed
        T = esc(:T)
        :($(_set(s)){$T})
    else
        _set(s)
    end
end
function _typedfun(s::SymbolFun)
    if s.typed
        T = esc(:T)
        :($(_fun(s)){$T})
    else
        _fun(s)
    end
end

# Base.lowercase is moved to Unicode.lowercase in Julia v0.7
using Unicode

_field(s::SymbolFS) = Symbol(replace(lowercase(string(s.s)), "." => "_"))

_getC(s::SymbolSet) = :(VectorOfConstraints{F,$(_typedset(s))})
_getC(s::SymbolFun) = _typedfun(s)

_callfield(f, s::SymbolFS) = :($f(model.$(_field(s))))
_broadcastfield(b, s::SymbolFS) = :($b(f, model.$(_field(s))))
_mapreduce_field(s::SymbolFS) = :(cur = op(cur, f(model.$(_field(s)))))
_mapreduce_constraints(s::SymbolFS) = :(cur = op(cur, f(model.$(_field(s)))))

"""
    struct_of_constraint_code(struct_name, types, field_types = nothing)

Given a vector of `n` `SymbolFun` or `SymbolSet` in `types`, creates a
struct of name `struct_name` that is a subtype of
`StructOfConstraint{T, C1, C2, ..., Cn}` if `field_types` is `nothing` and
a subtype of `StructOfConstraint{T}` otherwise.
It contains `n` field where the `i`th field has type `Ci` if `field_types` is
`nothing` and type `field_types[i]` otherwise.
If `types` is vector of `SymbolFun` (resp. `SymbolSet`) then the constraints
of that function (resp. set) type are stored in the corresponding field.
"""
function struct_of_constraint_code(struct_name, types, field_types = nothing)
    esc_struct_name = struct_name
    T = esc(:T)
    typed_struct = :($(esc_struct_name){$T})
    type_parametrized = field_types === nothing
    if type_parametrized
        field_types = [Symbol("C$i") for i in eachindex(types)]
        append!(typed_struct.args, field_types)
    end

    struct_def = :(struct $typed_struct <: StructOfConstraints end)

    for (t, field_type) in zip(types, field_types)
        field = _field(t)
        push!(struct_def.args[3].args, :($field::$field_type))
    end
    code = quote
        function $MOIU.broadcastcall(f::Function, model::$esc_struct_name)
            return $(Expr(:block, _callfield.(Ref(:f), types)...))
        end
        function $MOIU.broadcastvcat(f::Function, model::$esc_struct_name)
            return vcat($(_callfield.(Ref(:f), types)...))
        end
        function $MOIU.mapreduce_constraints(
            f::Function,
            op::Function,
            model::$esc_struct_name,
            cur,
        )
            return $(Expr(:block, _mapreduce_field.(types)...))
        end
    end

    for t in types
        if t isa SymbolFun
            fun = _fun(t)
            set = :(MOI.AbstractSet)
        else
            fun = :(MOI.AbstractFunction)
            set = _set(t)
        end
        field = _field(t)
        code = quote
            $code
            function $MOIU.constraints(
                model::$esc_struct_name,
                ::Type{<:$fun},
                ::Type{<:$set},
            ) where {S}
                return model.$field
            end
        end
    end
    supports_code = if eltype(types) <: SymbolFun
        quote
            function $MOI.supports_constraint(
                model::$esc_struct_name{$T},
                ::Type{F},
                ::Type{S},
            ) where {
                $T,
                F<:Union{$(_typedfun.(types)...)},
                S<:MOI.AbstractSet,
            }
                return $MOI.supports_constraint(constraints(model, F, S), F, S)
            end
        end
    else
        @assert eltype(types) <: SymbolSet
        quote
            function $MOI.supports_constraint(
                model::$esc_struct_name{$T},
                ::Type{F},
                ::Type{S},
            ) where {
                $T,
                F<:MOI.AbstractFunction,
                S<:Union{$(_typedset.(types)...)},
            }
                return $MOI.supports_constraint(constraints(model, F, S), F, S)
            end
        end
    end
    expr = Expr(:block, struct_def, supports_code, code)
    if !isempty(field_types)
        constructors = [:($field_type()) for field_type in field_types]
        # If there is no field type, the default constructor is sufficient and
        # adding this constructor will make a `StackOverflow`.
        constructor_code = :(function $typed_struct() where {$T}
            return $typed_struct($(constructors...))
        end)
        if type_parametrized
            append!(constructor_code.args[1].args, field_types)
        end
        push!(expr.args, constructor_code)
    end
    return expr
end


text:
using Test

include("dnd-character.jl")

@testset "Randomly generated ability is within range" begin
    for i in 1:1000
        @test 3 <= ability() <= 18
    end
end

@testset "Ability modifiers are integers" begin
    @test typeof(modifier(7)) <: Integer
end

@testset "Ability modifiers" begin
    @test modifier(3) == -4
    @test modifier(4) == -3
    @test modifier(5) == -3
    @test modifier(6) == -2
    @test modifier(7) == -2
    @test modifier(8) == -1
    @test modifier(9) == -1
    @test modifier(10) == 0
    @test modifier(11) == 0
    @test modifier(12) == 1
    @test modifier(13) == 1
    @test modifier(14) == 2
    @test modifier(15) == 2
    @test modifier(16) == 3
    @test modifier(17) == 3
    @test modifier(18) == 4
end

@testset "Randomly generated character is valid" begin
    # Helper method to check if all abilities of a character are within the expected ranges.
    function ischaracter(c)
        3 <= c.strength     <= 18 &&
        3 <= c.dexterity    <= 18 &&
        3 <= c.constitution <= 18 &&
        3 <= c.intelligence <= 18 &&
        3 <= c.wisdom       <= 18 &&
        3 <= c.charisma     <= 18 &&
        c.hitpoints == 10 + modifier(c.constitution)
    end

    for i in 1:1000
        c = DNDCharacter()
        @test ischaracter(c)
    end
end


text:
#!/usr/bin/env julia

# This is an example invocation of `debootstrap` to generate a Debian/Ubuntu-based rootfs
using Scratch, Tar, Pkg, Pkg.Artifacts, ghr_jll, SHA

if Sys.which("debootstrap") === nothing
    error("Must install `debootstrap`!")
end

# Utility functions
getuid() = ccall(:getuid, Cint, ())
getgid() = ccall(:getgid, Cint, ())

artifact_hash = create_artifact() do rootfs
    release = "buster"
    @info("Running debootstrap")
    run(`sudo debootstrap --variant=minbase --include=locales $(release) "$(rootfs)"`)

    # Remove special `dev` files
    @info("Cleaning up `/dev`")
    for f in readdir(joinpath(rootfs, "dev"); join=true)
        # Keep the symlinks around (such as `/dev/fd`), as they're useful
        if !islink(f)
            run(`sudo rm -rf "$(f)"`)
        end
    end

    # take ownership of the entire rootfs
    @info("Chown'ing rootfs")
    run(`sudo chown $(getuid()):$(getgid()) -R "$(rootfs)"`)

    # Write out a reasonable default resolv.conf
    open(joinpath(rootfs, "etc", "resolv.conf"), write=true) do io
        write(io, """
        nameserver 1.1.1.1
        nameserver 8.8.8.8
        nameserver 8.8.4.4
        nameserver 4.4.4.4
        """)
    end

    # Remove `_apt` user so that `apt` doesn't try to `setgroups()`
    @info("Removing `_apt` user")
    open(joinpath(rootfs, "etc", "passwd"), write=true, read=true) do io
        filtered_lines = filter(l -> !startswith(l, "_apt:"), readlines(io))
        truncate(io, 0)
        seek(io, 0)
        for l in filtered_lines
            println(io, l)
        end
    end

    # Set up the one true locale
    @info("Setting up UTF-8 locale")
    open(joinpath(rootfs, "etc", "locale.gen"), "a") do io
        println(io, "en_US.UTF-8 UTF-8")
    end
    @info("Regenerating locale")
    run(`sudo chroot --userspec=$(getuid()):$(getgid()) $(rootfs) locale-gen`)
    @info("Done!")
end

# Archive it into a `.tar.gz` file
@info("Archiving")
tarball_path = joinpath(@get_scratch!("archived"), "debian_minimal.tar.gz")
archive_artifact(artifact_hash, tarball_path)

# Hash the tarball
@info("Hashing tarball")
tarball_hash = open(io -> bytes2hex(sha256(io)), tarball_path)

# Upload it to `staticfloat/Sandbox.jl`
tag_name = "debian-minimal-$(bytes2hex(artifact_hash.bytes[end-3:end]))"
@info("Uploading to staticfloat/Sandbox.jl@$(tag_name)")
run(`$(ghr_jll.ghr()) -replace $(tag_name) $(tarball_path)`)

# Bind this artifact into our Artifacts.toml
tarball_url = "https://github.com/staticfloat/Sandbox.jl/releases/download/$(tag_name)/$(basename(tarball_path))"
bind_artifact!(
    joinpath(dirname(@__DIR__), "Artifacts.toml"),
    "debian-minimal-rootfs",
    artifact_hash;
    download_info=[(tarball_url, tarball_hash)],
    lazy=true,
    force=true,
)


text:
#!/usr/bin/env julia

using DelimitedFiles

using WaterMeanMeanForce

using ArgParse

s = ArgParseSettings()
s.autofix_names = true

models = join(list_models(), ", ")
integrators = join(list_integrators(), ", ")

@add_arg_table! s begin
    "--model"
        metavar = "M"
        range_tester = in(list_models())
        help = "water model ($(models))"
        required = true
    "--integrator"
        metavar = "I"
        range_tester = in(list_integrators())
        help = "molecular dynamics integrator ($(integrators))"
        required = true
    "-T"
        metavar = "T"
        help = "temperature (K)"
        arg_type = Float64
        required = true
    "-R"
        metavar = "R"
        help = "constraint/restraint distance (nm)"
        arg_type = Float64
        required = true
    "--num-links"
        metavar = "P"
        help = "number of Trotter links"
        arg_type = Int
        required = true
    "--time-step"
        metavar = "T"
        help = "time step (ps)"
        arg_type = Float64
        required = true
    "--centroid-friction"
        metavar = "G"
        help = "centroid friction (1 / ps)"
        arg_type = Float64
        required = true
    "--equil-duration"
        metavar = "D"
        help = "duration of equilibration (ps)"
        arg_type = Float64
        required = true
    "--umbrella-k"
        metavar = "K"
        help = "umbrella sampling force constant (kJ / nm^2 mol)"
        arg_type = Float64
    "--initial-configuration"
        metavar = "FILE"
        help = "path to initial configuration input file"
    "--final-configuration"
        metavar = "FILE"
        help = "path to final configuration output file"
end

c = parse_args(ARGS, s, as_symbols=true)

if !isnothing(c[:initial_configuration])
    qs_init = dropdims(readdlm(c[:initial_configuration]); dims=2)
else
    qs_init = nothing
end

mdp_equil = MolecularDynamicsParameters(c[:centroid_friction], c[:time_step],
                                        c[:equil_duration])

if isnothing(c[:umbrella_k])
    s = run_constrained_equilibration(get_model(c[:model]),
                                      get_integrator(c[:integrator]), 2, c[:T],
                                      c[:R], c[:num_links], mdp_equil; qs_init)
else
    s = run_restrained_equilibration(get_model(c[:model]),
                                     get_integrator(c[:integrator]), 2, c[:T],
                                     c[:R], c[:num_links], mdp_equil,
                                     c[:umbrella_k]; qs_init)
end

if !isnothing(c[:final_configuration])
    writedlm(c[:final_configuration], s.qs)
end


text:
#
# FourierOptics.jl -
#
# Tools for Fourier optics computations in Julia.
#

__precompile__(true)

module FourierOptics

export
    CoordinateTransform,
    Region,
    center,
    circularmask!,
    circularmask,
    fftshiftphasor!,
    fftshiftphasor,
    recenter

using LazyAlgebra
import LazyAlgebra: Complexes, Reals

using AbstractFFTs
import FFTW

# Declaration of functions to be extended by sub-modules.
function center end
function recenter end

include("units.jl")
using .Units

include("coords.jl")
import .CoordinateTransforms:
    CoordinateTransform,
    compose,
    jacobian,
    origin

include("regions.jl")
using .Regions

include("types.jl")
include("utils.jl")
include("lenses.jl")

end # module


text:
using SteadyStateDiffEq, DiffEqBase, NLsolve, Sundials
using Test

function f(du,u,p,t)
  du[1] = 2 - 2u[1]
  du[2] = u[1] - 4u[2]
end
u0 = zeros(2)
prob = SteadyStateProblem(f,u0)
abstol = 1e-8
sol = solve(prob,SSRootfind())
@test sol.retcode == :Success
p = nothing

du = zeros(2)
f(du,sol.u,nothing,0)
@test maximum(du) < 1e-11

prob = ODEProblem(f,u0,(0.0,1.0))
prob = SteadyStateProblem(prob)
sol = solve(prob,SSRootfind(nlsolve = (f,u0,abstol) -> (res=NLsolve.nlsolve(f,u0,autodiff=:forward,method=:newton,iterations=Int(1e6),ftol=abstol);res.zero) ))
@test sol.retcode == :Success

f(du,sol.u,nothing,0)
@test du == [0,0]

# Use Sundials
sol = solve(prob,SSRootfind(nlsolve = (f,u0,abstol) -> (res=Sundials.kinsol(f,u0)) ))
@test sol.retcode == :Success
f(du,sol.u,nothing,0)
@test du == [0,0]

using OrdinaryDiffEq
sol = solve(prob,DynamicSS(Rodas5()))
@test sol.retcode == :Success

f(du,sol.u,p,0)
@test du ≈ [0,0] atol = 1e-7

sol = solve(prob,DynamicSS(Rodas5(),tspan=1e-3))
@test sol.retcode != :Success

sol = solve(prob,DynamicSS(CVODE_BDF()),dt=1.0)
@test sol.retcode == :Success

f(du,sol.u,p,0)
@test du ≈ [0,0] atol = 1e-6


text:
@testset "2d approximation" begin
  l = collect(-2:0.2:2);
  ins = transpose(hcat(collect.(collect(product(l, l))[:])...));
  func = (x, y) -> sin(x) + exp(-y^2)
  fs = [func(ins[i, 1], ins[i, 2]) for i in 1:size(ins, 1)]

  m = mwls_kd(ins, fs, 2, (d, eps) -> exp(-d^2), maxDegree = 3)
  cll = mwls_cll(hcat(ins, fs), 2, (d, eps) -> exp(-d^2), maxDegree = 3)
  teststatus = true
  for i in size(ins, 1)
    teststatus &= isapprox(m([ins[i, 1], ins[i, 2]]),
                           func(ins[i, 1], ins[i, 2]),
                           atol = 0.2)
  end
  @test teststatus == true

  for i in size(ins, 1)
    teststatus &= isapprox(cll([ins[i, 1], ins[i, 2]]),
                           func(ins[i, 1], ins[i, 2]),
                           atol = 0.3)
  end
  @test teststatus == true
end


text:
module Darwin

export AbstractEvolutionaryModel,
    FitnessValue

abstract type AbstractPopulationModel{G} end

const FitnessValue = Float64

include("utils.jl")
include("individual.jl")
include("rates.jl")
include("fitness.jl")
include("population_model.jl")

include("selection.jl")
include("mutation.jl")
include("crossover.jl")

include("GA.jl")

end # module


text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder

# Collection of sources required to build SymEngine
name = "SymEngine"
version = v"0.6.0"
sources = [
    ArchiveSource("https://github.com/symengine/symengine/releases/download/v$(version)/symengine-$(version).tar.gz",
                  "4d2caa86c03eaaa8ed004084d02f87b5c51b6229f8ba70d161227e22d6302f0a"),
]

# Bash recipe for building across all platforms

script = raw"""
cd $WORKSPACE/srcdir/symengine-*
mkdir build
cd build
cmake -DCMAKE_INSTALL_PREFIX=${prefix} \
      -DCMAKE_TOOLCHAIN_FILE=${CMAKE_TARGET_TOOLCHAIN} \
      -DBUILD_TESTS=no \
      -DBUILD_BENCHMARKS=no \
      -DBUILD_SHARED_LIBS=yes \
      -DWITH_MPC=yes \
      -DWITH_COTIRE=no \
      -DWITH_SYMENGINE_THREAD_SAFE=yes ..
make -j${nproc}
make install
"""

platforms = expand_cxxstring_abis(supported_platforms())

# The products that we will ensure are always built
products = [
    LibraryProduct(["libsymengine", "libsymengine-$(version.major).$(version.minor)"], :libsymengine)
]

# Dependencies that must be installed before this package can be built
dependencies = [
    Dependency("GMP_jll", v"6.1.2"),
    Dependency("MPFR_jll", v"4.0.2"),
    "MPC_jll",
]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies)

max_stars_repo_name: jameskermode/Yggdrasil

max_issues_repo_name: jameskermode/Yggdrasil

max_forks_repo_name: jameskermode/Yggdrasil


text:
function qualityMeasure(xN::Array{Int64,1}, yN::Array{Int64,1},
                        xU::Array{Int64,1}, yU::Array{Int64,1}
                       )

    # ---- Set the reference point
    xRef = xN[end] + 1
    yRef = yN[1] + 1

    # ---- Measure the area corresponding to the set of non-dominated points Y_N
    areaN = 0
    for i = 1:length(xN)-1
        areaN = areaN + (xN[i+1]-xN[i]) * (yRef - yN[i])
    end
    areaN = areaN + (xRef-xN[end]) * (yRef - yN[end])

    # ---- Measure the (restricted) area corresponding to the bound set U
    areaU = 0
    for i = 1:length(xU)-1
        areaU = areaU + ( min(xRef,xU[i+1]) - min(xRef,xU[i]) ) * (yRef - min(yRef,yU[i]) )
    end
    areaU = areaU + ( xRef - min(xRef,xU[end]) ) * (yRef - min(yRef,yU[end]) )

    # ---- Computer the ratio
    result = areaU/areaN
    return result

end

max_stars_repo_name: ducanh-le/GravityMachine

max_issues_repo_name: ducanh-le/GravityMachine

max_forks_repo_name: ducanh-le/GravityMachine


text:
using Roots

import Base: show

## Table 1 from TOMS748 by Alefeld, Potra, Shi

type Func
    name :: Symbol
    val :: Function
    bracket :: Function
    params :: Vector{Any}
end
function show(io::IO, f::Func)
    @printf io "Func(%s)" f.name
end
known_functions = Func[]

## Construct a function object, and check root brackets
macro Func(name)
    quote
        f = Func($name, val, bracket, params)
        for p in params
            b = bracket(p)
            @assert val(p, b[1]) * f.val(p, b[2]) < 0 "Invalid bracket"
        end
        push!(known_functions, f)
        f
    end
end

func1 = let
    val(_, x) = sin(x) - x/2
    bracket(_) = [0.5pi, pi]
    params = [()]
    @Func :func1
end

func2 = let
    val(n, x) = -2*sum([(2i-5)^2/(x-i*i)^3 for i=1:20])
    bracket(n) = [n^2+1e-9, (n+1)^2-1e-9]
    params = [1:10]
    @Func :func2
end

func3 = let
    val(p, x) = p[1]*x*exp(p[2]*x)
    bracket(p) = [-9., 31.]
    params = [(-40.,-1.), (-100., -2.), (-200., -3.)]
    @Func :func3
end

func4 = let
    val(p, x) = x^p[2] - p[1]
    bracket(p) = p[3]
    params = (Float64, Float64, Vector{Float64})[]
    for a in [0.2, 1.], n in 4:2:12
        push!(params, (a, n, [0., 5.]))
    end
    for n in 8:2:14
        push!(params, (1., n, [-0.95, 4.05]))
    end
    @Func :func4
end

func5 = let
    val(p, x) = sin(x) - 0.5
    bracket(p) = [0., 1.5]
    params = [()]
    @Func :func5
end

func6 = let
    val(n, x) = 2x*exp(-n)-2exp(-n*x)+1.
    bracket(n) = [0., 1.]
    params = vcat(1:5, 20:20:100)
    @Func :func6
end

func7 = let
    val(n, x) = (1+(1-n)^2)*x-(1-n*x)^2
    bracket(n)= [0., 1.]
    params = [5., 10., 20.]
    @Func :func7
end

func8 = let
    val(n, x) = x^2-(1-x)^n
    bracket(n) = [0., 1.]
    params = [2., 5., 10., 15., 20.]
    @Func :func8
end

func9 = let
    val(n, x) = (1+(1-n)^4)*x-(1-n*x)^4
    bracket(n) = [0., 1.]
    params = [1., 2., 4., 5., 8., 15., 20.]
    @Func :func9
end

func10 = let
    val(n, x) = exp(-n*x)*(x-1) + x^n
    bracket(n) = [0., 1.]
    params = [1, 5, 10, 15, 20]
    @Func :func10
end

func11 = let
    val(n, x) = (n*x-1)/((n-1)*x)
    bracket(n) = [0.01, 1.]
    params = [2, 5, 15, 20]
    @Func :func11
end

func12 = let
    val(n, x) = x^(1/n)-n^(1/n)
    bracket(n) = [1., 100.]
    params = vcat(2:6, 7:2:33)
    @Func :func12
end

func13 = let
    val(n, x) = x == 0. ? 0. : x/exp(1/(x*x))
    bracket(n) = [-1., 4.]
    params = [()]
    @Func :func13
end

func14 = let
    val(n, x) = x >= 0 ? n/20*(x/1.5+sin(x)-1) : -n/20
    bracket(n) = [-1e4, 0.5pi]
    params = [1:40]
    @Func :func14
end

func15 = let
    function val(n, x)
        if x > 2e-3/(1+n)
            e - 1.859
        elseif x < 0
            -0.859
        else
            exp(0.5e3(n+1)x)-1.859
        end
    end
    bracket(n) = [-1e4, 1e-4]
    params = vcat(20:40, 100:100:1000)
    @Func :func15
end

type MethodResults
    name
    evalcount :: Int
    maxresidual :: Float64
    failures :: Vector{(Func, Int)}
end
MethodResults() = MethodResults(nothing, 0, 0., (Func, Int)[])
show(io::IO, results::MethodResults) =
    print(io, "MethodResults($(results.name), evalcount=$(results.evalcount), numfailures=$(length(results.failures)), maxresidual=$(results.maxresidual))")

## Run a method on all known functions.
function run_tests(method; verbose=false, trace=false, name=nothing, abandon=false)
    results = MethodResults()
    results.name = name
    for f in known_functions
        for i in 1:length(f.params)
            p = f.params[i]
            evalcount = 0
            function feval(x)
                evalcount += 1
                result = f.val(p, x)
                trace && @printf "%s[%d]: %s ⇒ %s\n" f i x result
                result
            end
            result, residual = nothing, nothing
            try
                result = method(feval, f.bracket(p))
                residual = f.val(p, result)
                verbose && @printf "%s[%d] ⇒ %d / %s, residual %.5e\n" f i evalcount result residual
            catch ex
                verbose && @printf "%s[%d] ⇒ FAILED: %s\n" f i ex
                push!(results.failures, (f, i))
                abandon && rethrow(ex)
            end
            results.evalcount += evalcount
            ## Some functions might return non-real values on failures
            if isa(result, Float64) && isa(residual, Float64) && isfinite(residual)
                results.maxresidual = max(results.maxresidual, abs(residual))
            end
        end
    end
    results
end

function run_tests()
    @printf "%s\n" run_tests(fzero, name="fzero")
    @printf "%s\n" run_tests((f, b) -> fzero(f, mean(b)), name="fzero(x0)")
    @printf "%s\n" run_tests((f, b) -> fzero(f, b), name="fzero(b)")
    @printf "%s\n" run_tests((f, b) -> fzero(f, mean(b), b), name="fzero(x0,b)")
    @printf "%s\n" run_tests((f, b) -> fzero(f, mean(b), order=1), name="fzero/1")
    @printf "%s\n" run_tests((f, b) -> fzero(f, mean(b), order=2), name="fzero/2")
    @printf "%s\n" run_tests((f, b) -> fzero(f, mean(b), order=5), name="fzero/5")
    @printf "%s\n" run_tests((f, b) -> fzero(f, mean(b), order=8), name="fzero/8")
    @printf "%s\n" run_tests((f, b) -> fzero(f, mean(b), order=16), name="fzero/16")
    @printf "%s\n" run_tests((f, b) -> newton(f, mean(b)), name="newton")
    @printf "%s\n" run_tests((f, b) -> halley(f, mean(b)), name="halley")
    @printf "%s\n" run_tests((f, b) -> secant_method(f, b[1], b[2]), name="secant")
end

nothing


text:
@inline function Vreman_eddy_viscosity(S::SymTen,ω::Vec,c::Real,Δ2::Real)
    α = S + AntiSymTen(-0.5*ω)
    β = Δ2 * symmetric(α' ⋅ α)
    B = β.xx*(β.yy + β.zz) + β.yy*β.zz - (β.xy^2 + β.xz^2 + β.yz^2)
    B = ifelse(B<0.0,0.0,B)
    return c*fsqrt(B/(α:α))
end


text:
import Base: Pair, Tuple, show, ==
import LightGraphs: AbstractEdge, src, dst, reverse

abstract type AbstractSimpleWeightedEdge{T} <: AbstractEdge{T} end

struct SimpleWeightedEdge{T<:Integer, U<:Real} <: AbstractSimpleWeightedEdge{T}
    src::T
    dst::T
    weight::U
end

SimpleWeightedEdge(t::NTuple{2}) = SimpleWeightedEdge(t[1], t[2], one(Float64))
SimpleWeightedEdge(t::NTuple{3}) = SimpleWeightedEdge(t[1], t[2], t[3])
SimpleWeightedEdge(p::Pair) = SimpleWeightedEdge(p.first, p.second, one(Float64))
SimpleWeightedEdge{T, U}(p::Pair) where T<:Integer where U <: Real = SimpleWeightedEdge(T(p.first), T(p.second), one(U))
SimpleWeightedEdge{T, U}(t::NTuple{3}) where T<:Integer where U <: Real = SimpleWeightedEdge(T(t[1]), T(t[2]), U(t[3]))
SimpleWeightedEdge{T, U}(t::NTuple{2}) where T<:Integer where U <: Real = SimpleWeightedEdge(T(t[1]), T(t[2]), one(U))
SimpleWeightedEdge(x, y) = SimpleWeightedEdge(x, y, one(Float64))
eltype(e::T) where T<:AbstractSimpleWeightedEdge= eltype(src(e))

# Accessors
src(e::AbstractSimpleWeightedEdge) = e.src
dst(e::AbstractSimpleWeightedEdge) = e.dst
weight(e::AbstractSimpleWeightedEdge) = e.weight

# I/O
show(io::IO, e::AbstractSimpleWeightedEdge) = print(io, "Edge $(e.src) => $(e.dst) with weight $(e.weight)")

# Conversions
Tuple(e::AbstractSimpleWeightedEdge) = (src(e), dst(e), weight(e))

(::Type{SimpleWeightedEdge{T, U}})(e::AbstractSimpleWeightedEdge) where {T <: Integer, U <: Real} = SimpleWeightedEdge{T, U}(T(e.src), T(e.dst), U(e.weight))

# Convenience functions - note that these do not use weight.
reverse(e::T) where T<:AbstractSimpleWeightedEdge = T(dst(e), src(e), weight(e))
==(e1::AbstractSimpleWeightedEdge, e2::AbstractSimpleWeightedEdge) = (src(e1) == src(e2) && dst(e1) == dst(e2))


text:
# TODO: maybe more fancy file format and correctness checking should be done

"""
    hg_save(io::IO, h::Hypergraph)

Saves a hypergraph `h` to an output stream `io`.

"""
function hg_save(io::IO, h::Hypergraph)
    println(io, length(h.v2he), " ", length(h.he2v))
    for he in h.he2v
        skeys = sort(collect(keys(he)))
        println(io, join(["$k=$(he[k])" for k in skeys], ' '))
    end
end

"""
    hg_save(fname::AbstractString, h::Hypergraph)

Saves a hypergraph `h` to a file `fname`.

"""
hg_save(fname::AbstractString, h::Hypergraph) =
    open(io -> hg_save(io, h), fname, "w")

"""
    hg_load(fname::AbstractString, T::Type{<:Real})

Loads a hypergraph from a stream `io`. The second argument
`T` represents type of data in the hypegraph.

Skips an initial comment.

"""
function hg_load(io::IO, T::Type{<:Real})
    line = readline(io)

    if startswith(line, "\"\"\"")
      singleline = true
        while(
            !( (!singleline && endswith(line, "\"\"\"")) ||
            (singleline && endswith(line, "\"\"\"") && length(line)>5)
            ) &&
            !eof(io)
            )
                line = readline(io)
                singleline = false
        end
        if eof(io)
            throw(ArgumentError("malformed input"))
        end
       line = readline(io)
    end

    l = split(line)
    length(l) == 2 || throw(ArgumentError("expected two integers"))
    n, k = parse.(Int, l)
    h = Hypergraph{T}(n, k)
    lastv = 0
    for i in 1:k
        for pos in split(readline(io))
            entry = split(pos, '=')
            length(entry) == 2 || throw(ArgumentError("expected vertex=weight"))
            v = parse(Int, entry[1])
            w = parse(T, entry[2])
            if v > lastv
                lasti = v
            else
                throw(ArgumentError("vertices in hyperedge must be sorted"))
            end
            h[v, i] = w
        end
    end
    # we ignore lines beyond k+1 in the file
    h
end

"""
    hg_load(fname::AbstractString, T::Type{<:Real})

Loads a hypergraph from a file `fname`. The second argument
`T` represents type of data in the hypegraph

"""
hg_load(fname::AbstractString, T::Type{<:Real}) =
    open(io -> hg_load(io, T), fname, "r")


text:
#Helper function used due to performance bug in @threads.
function _loopbody!(
    pivot::U,
    nvg::U,
    dists::Matrix{T},
    parents::Matrix{U}
    ) where T<:Real where U<:Integer
    # Relax dists[u, v] = min(dists[u, v], dists[u, pivot]+dists[pivot, v]) for all u, v
    @inbounds @threads for v in one(U):nvg
        d = dists[pivot, v]
        if d != typemax(T) && v != pivot
            p = parents[pivot, v]
            @inbounds for u in one(U):nvg
                ans = (dists[u, pivot] == typemax(T) || u == pivot ? typemax(T) : dists[u, pivot] + d)
                if dists[u, v] > ans
                    dists[u, v] = ans
                    parents[u, v] = p
                end
            end
        end
    end
end

function floyd_warshall_shortest_paths(
    g::AbstractGraph{U},
    distmx::AbstractMatrix{T}=weights(g)
) where T<:Real where U<:Integer
    nvg = nv(g)
    dists = fill(typemax(T), (Int(nvg), Int(nvg)))
    parents = zeros(U, (Int(nvg), Int(nvg)))

    for v in 1:nvg
        dists[v, v] = zero(T)
    end
    undirected = !is_directed(g)
    for e in edges(g)
        u = src(e)
        v = dst(e)

        d = distmx[u, v]

        dists[u, v] = min(d, dists[u, v])
        parents[u, v] = u
        if undirected
            dists[v, u] = min(d, dists[v, u])
            parents[v, u] = v
        end
    end

    for pivot in vertices(g)
        _loopbody!(pivot, nvg, dists, parents) #Due to bug in @threads
    end
    fws = FloydWarshallState(dists, parents)
    return fws
end


text:
using Test
using CliOptions

include("testutils.jl")

@testset "parse_args()" begin
    @testset "type of args; $(typeof(v))" for v in [
        ["-abc", "foo", "bar"],
        ("-abc", "foo", "bar"),
        split("-abc foo bar"),
    ]
        args = v
        spec = CliOptionSpec(FlagOption("-a"),
                             CounterOption("-b"),
                             Option("-c"),
                             Positional("d"),
                             onerror = error)
        result = parse_args(spec, args)
        @test result.a == true
        @test result.b == 1
        @test result.c == "foo"
        @test result.d == "bar"
    end

    @testset "normalization; $(v[1])" for v in [
        (["-fa", "foo"], [FlagOption("-f"), Option("-a")], [:f => true, :a => "foo"]),
        (["-af", "foo"], [FlagOption("-f"), Option("-a")], ErrorException),
        (["--foo-bar=baz"], [Option("--foo-bar")], [:foo_bar => "baz"]),
    ]
        args, options, expected = v
        spec = CliOptionSpec(options...; onerror = error)
        if expected isa Type
            @test_throws expected parse_args(spec, args)
        else
            args = parse_args(spec, args)
            for (k, v) in expected
                @test getproperty(args, k) == v
            end
        end
    end

    @testset "Option" begin
        @testset "required" begin
            spec = CliOptionSpec(
                Option("-a"),
                onerror = error,
            )
            @test_throws ErrorException parse_args(spec, String[])
            @test_throws ErrorException parse_args(spec, ["-a"])
            args = parse_args(spec, ["-a", "bar"])
            @test args.a == "bar"
            @test_throws ErrorException parse_args(spec, ["-a", "bar", "-a"])
        end

        @testset "omittable" begin
            spec = CliOptionSpec(
                Option("-a"; default = nothing),
                onerror = error,
            )
            args = parse_args(spec, String[])
            @test args.a === nothing
        end
    end

    @testset "FlagOption; $(v[1])" for v in [
        (["-a"], true),
        (["-b"], false),
        ([], false),
    ]
        args, expected = v
        spec = CliOptionSpec(
            FlagOption("-a"; negators = "-b"),
            onerror = error,
        )
        result = parse_args(spec, args)
        @test result.a == expected
        @test result.b == !expected
    end

    @testset "CounterOption; $(v[1:3])" for v in [
        (Int8, 127, "--add", Int8(-128)),  # overflow
        (Int8, 126, "--add", Int8(127)),
        (Int8, -127, "--sub", Int8(-128)),
        (Int8, -128, "--sub", Int8(127)),  # underflow
        (Int128, 0, "--add --sub", Int128(0)),
    ]
        T, default, args, expected = v
        spec = CliOptionSpec(CounterOption(T, "--add";
                                           decrementers = ["--sub"],
                                           default = default),
                             onerror = error)
        if expected isa Type
            @test_throws expected parse_args(spec, split(args))
        else
            args = parse_args(spec, split(args))
            @test propertynames(args) == [:add]
            @test args.add == expected
            @test typeof(args.add) == typeof(expected)
        end
    end

    @testset "HelpOption; $(v[1])" for v in [
        ("-h", (true, 1)),
        ("-a -h", (false, 0)),
        ("-a -b -h", (true, 1)),
    ]
        args, expected = v
        counter = 0
        spec = CliOptionSpec(
            Option("-a"; default = "foo"),
            HelpOption(),
            onhelp = () -> counter += 1,
            onerror = error,
        )
        result = parse_args(spec, split(args))
        @test result.help == expected[1]
        @test counter == expected[2]
    end

    @testset "HelpOption; onhelp = $(v[1])" for v in [
        ("Integer", 42, (c) -> error("foobar$c"), (ErrorException, "foobar42")),
        ("Nothing", nothing, Base.exit, (true, nothing)),
        ("Function", () -> error("foo"), (c) -> error("$c"), (ErrorException, "foo")),
    ]
        _, onhelp, exitfunc, expected = v
        spec = CliOptionSpec(
            HelpOption(),
            onhelp = onhelp,
            onerror = error,
        )
        CliOptions._mock_exit_function(exitfunc) do
            buf = IOBuffer()
            redirect_stdout(buf) do
                if expected[1] isa Type
                    @test_throws expected[1] parse_args(spec, ["-h"])
                    output = String(take!(buf))
                    @test occursin(output, expected[2])
                else
                    result = parse_args(spec, ["-h"])
                    @test result.help == expected[1]
                end
            end
        end
    end

    @testset "Positional; $(v[1]), $(v[4])" for v in [
        # single, required
        ("single, required", false, missing, String[], ErrorException),
        ("single, required", false, missing, ["a"], "a"),
        ("single, required", false, missing, ["a", "b"], ErrorException),
        ("single, required", false, missing, ["-1"], "-1"),
        ("single, required", false, missing, ["-7"], ErrorException),
        ("single, required", false, missing, ["-a"], ErrorException),

        # single, omittable
        ("single, omittable", false, "foo.txt", String[], "foo.txt"),
        ("single, omittable", false, "foo.txt", ["a"], "a"),
        ("single, omittable", false, "foo.txt", ["a", "b"], ErrorException),

        # multiple, required
        ("multiple, required", true, missing, String[], ErrorException),
        ("multiple, required", true, missing, ["a"], ["a"]),
        ("multiple, required", true, missing, ["a", "b"], ["a", "b"]),
        ("multiple, required", true, missing, ["a", "-1"], ["a", "-1"]),
        ("multiple, required", true, missing, ["a", "-7"], ["a"]),
        ("multiple, required", true, missing, ["a", "-a"], ErrorException),
        ("multiple, required", true, missing, ["-1", "a"], ["-1", "a"]),
        ("multiple, required", true, missing, ["-7", "a"], ["a"]),
        ("multiple, required", true, missing, ["-a", "a"], ErrorException),

        # multiple, omittable
        ("multiple, omittable", true, "foo.txt", String[], "foo.txt"),
        ("multiple, omittable", true, "foo.txt", ["a"], ["a"]),
        ("multiple, omittable", true, "foo.txt", ["a", "b"], ["a", "b"]),
    ]
        title, multiple, default, args, expected = v
        spec = CliOptionSpec(
            FlagOption("-7"),
            Positional("file", "files"; multiple = multiple, default = default),
            onerror = error,
        )
        if expected == ErrorException
            @test_throws expected parse_args(spec, args)
        else
            result = parse_args(spec, args)
            @test result.file == expected
            @test result.files == expected
        end
    end

    @testset "OptionGroup" begin
        spec = CliOptionSpec(
            OptionGroup(
                Option("-a"),
                Option("-b"),
            ),
            onerror = error,
        )
        @test_throws ErrorException parse_args(spec, String[])
        @test_throws ErrorException parse_args(spec, split("-a foo"))
        @test_throws ErrorException parse_args(spec, split("-b bar"))
        args = parse_args(spec, split("-a foo -b bar"))
        @test args.a == "foo"
        @test args.b == "bar"
        @test_throws ErrorException parse_args(spec, split("-a foo -b bar baz"))
    end

    @testset "MutexGroup" begin
        spec = CliOptionSpec(
            MutexGroup(
                Option("-a"),
                Option("-b"),
                Option("-c"),
            ),
            onerror = error,
        )
        @test_throws ErrorException parse_args(spec, String[])
        args = parse_args(spec, split("-a foo"))
        @test args.a == "foo"
        @test args.b === missing
        @test args.c === missing
        args = parse_args(spec, split("-b bar"))
        @test args.a === missing
        @test args.b == "bar"
        @test args.c === missing
        @test_throws ErrorException parse_args(spec, split("-a foo -b bar"))
        @test_throws ErrorException parse_args(spec, split("-a foo -c baz"))
        @test_throws ErrorException parse_args(spec, split("-a foo -b bar -c baz"))
        @test_throws ErrorException parse_args(spec, split("-a foo quux"))
    end

    @testset "onerror = $(v[1])" for v in [
        ("Integer", 42, (c) -> error("foobar$c"), (ErrorException, "foobar42")),
        ("Nothing", nothing, () -> error("foo bar"), ("\"-a\"", nothing)),
        ("Function", (msg) -> error("foo"), (c) -> error("$c"), (ErrorException, "foo")),
    ]
        t, onerror, exitfunc, expected = v
        spec = CliOptionSpec(
            FlagOption("-f"),
            Option("-a"),
            CounterOption("-v"),
            onerror = onerror
        )
        CliOptions._mock_exit_function(exitfunc) do
            buf = IOBuffer()
            redirect_stderr(buf) do
                if expected[1] isa Type
                    tr = @test_throws expected[1] parse_args(spec, ["-fva"])
                    if tr isa Test.Pass
                        @test occursin(expected[2], tr.value.msg)
                    end
                else
                    result = parse_args(spec, ["-fva"])
                    @test occursin(expected[1], result._errors[1])
                end
            end
        end
    end

    @testset "use_double_dash; $(v[1]), $(v[2])" for v in [
        (false, "-a foo -- --", ["foo", "--", "--"]),
        (false, "-a -- foo --", ["--", "foo", "--"]),
        (false, "-- -a foo --", ["foo", "--"]),  # stops at `-a`, overwrites ["--"] with ["foo"]
        (true, "-a foo -- --", ["foo", "--"]),
        (true, "-a -- foo --", ["foo", "--"]),
        (true, "-- -a foo --", ["-a", "foo", "--"]),
    ]
        use_double_dash, args, expected = v
        spec = CliOptionSpec(
            FlagOption("-a"),
            Positional("file", "files"; multiple = true),
            use_double_dash = use_double_dash,
            onerror = error,
        )
        result = parse_args(spec, split(args))
        @test result.files == expected
    end

    @testset "ParseResult._errors, $(v[1])" for v in [
        ("1", split("-n 1 -x"), [
            "Unrecognized argument: \"-x\"",
        ]),
        ("2", split("-n a -x"), [
            "Invalid value for -n: \"a\" (must be one of \"1\", \"2\" or \"3\")",
            "Unrecognized argument: \"a\"",
            "Unrecognized argument: \"-x\"",
        ]),
    ]
        t, args, expected = v
        spec = CliOptionSpec(
            Option("-n"; requirement = ["1", "2", "3"]),
            onerror = nothing,
        )
        result = parse_args(spec, args)
        @test result._errors == expected
    end
end


text:
const k_pio2 = Double64(1.5707963267948966, 6.123233995736766e-17)

function atan(x::DoubleFloat{T}) where {T<:IEEEFloat}
   signbit(x) && return -atan(abs(x))
   isinf(x) && return DoubleFloat{T}(pi)/2
   iszero(x) && return x
   y = DoubleFloat{T}(atan(x.hi))
   s, c = sin(y), cos(y)
   c2 = square(c)
   t = s/c
   z = y - (t - x)*c2 # z - (tan(x)-x)*(cos(x)^2)
   s, c = sin(z), cos(z)
   t = s/c
   z = z - (t - x)*c2 # z - (tan(x)-x)*(cos(x)^2)
   return z
end

function atan(y::DoubleFloat{T}, x::DoubleFloat{T}) where {T<:IEEEFloat}
   iszero(y) && iszero(x) && return x
   iszero(x) && return copysign(pio2(DoubleFloat{T}), y)
   iszero(y) && return signbit(x) ? pio1(DoubleFloat{T}) : zero(DoubleFloat{T})
   atanyx = atan(y/x)
   if x > 0
      atanyx
   elseif x < 0
      if y >= 0
         atanyx + T(pi)
      else
         atanyx - T(pi)
      end
   elseif y > 0
       T(pi)/2
   elseif y < 0
      -T(pi)/2
   else
      zero(T)
   end
end

function asin(x::DoubleFloat{T}) where {T<:IEEEFloat}
   abs(x) > 1.0 && throw(DomainError("$x"))
   signbit(x) && return -asin(abs(x))
   y = x
   y = y / (1.0 + sqrt(1.0 - square(y)))
   z = atan(y)
   return DoubleFloat{T}(z.hi+z.hi, z.lo+z.lo)
end

function acos(x::DoubleFloat{T}) where {T<:IEEEFloat}
   abs(x) > 1.0 && throw(DomainError("$x"))
   signbit(x) && return DoubleFloat{T}(onepi - acos(abs(x)))
   y = x
   y = sqrt(1.0 - square(y)) / (1.0 + y)
   z = atan(y)
   return DoubleFloat{T}(z.hi+z.hi, z.lo+z.lo)
end

acsc(x::DoubleFloat{T}) where {T<:IEEEFloat} = asin(inv(x))
asec(x::DoubleFloat{T}) where {T<:IEEEFloat} = acos(inv(x))

function acot(x::DoubleFloat{T}) where {T<:IEEEFloat}
   signbit(x) && return DoubleFloat{T}(onepi - acot(abs(x)))
   iszero(x) && return k_pio2
   z = k_pio2 - atan(abs(x))
   if signbit(x.hi)
      z = -z
   end
   return z
end


text:
#' # Representing data in a computer
#' The core of data science and machine learning is **data**: we are interested in extracting knowledge from data.
#'
#' But how exactly do computers represent data? Let's find out exactly what an "artificial intelligence" has at its disposal to learn from.
import JuliaAcademyData
JAD = JuliaAcademyData

#' ## Data is represented as arrays
#'
#' Let's take a look at some fruit. Using the `Images.jl` library, we can load in some images:
import Images

#' Load apple
apple = Images.load("data/10_100.jpg")
#' Load banana
banana = Images.load("data/104_100.jpg")

#' Here we have images of apples and bananas. We would eventually like to build a program that can automatically distinguish between the two. However, the computer doesn't "see" an apple or a banana; instead, it just sees numbers.
#'
#' An image is encoded in something called an **array**, which is like a container that has boxes or slots for individual pieces of data:
#'
#' An array is a bunch of numbers in connected boxes; the figure above shows a 1-dimensional array. Our images are instead 2-dimensional arrays, or matrices, of numbers, arranged something like this:
#'
#'
#' ![array](data/array2d_s.png)
#'

typeof(apple)
#'
size(apple)
#' We can grab the datum stored in the box at row `i` and column `j` by *indexing* using square brackets: `[i, j]`. For example, let's get the pixel (piece of the image) in box $(40, 60)$, i.e. in the 40th row and 60th column of the image:
apple
#'
dump(apple[40,60])
#'
x = apple[40,60]
#'
apple[18:20,29:31]

#' We see that Julia displays a coloured box! Julia, via the `Colors.jl` package, is clever enough t'o display colours in a way that is useful to us humans!
#'
#' So, in fact, an image is a 2D array, in which each element of the array is an object (a collection of numbers) describing a coloured pixel.

#' ## Colors as numbers
#'
#' How, then, are these colors actually stored? Computers store colors in RGB format, that is they s'tore a value between 0 and 1 for each of three "channels": red, green, and blue. Here, 0 means n'one of that color and 1 means the brightest form of that color. The overall color is a c'ombination of those three colors.
#'
#' For example, we can pull out the `red` value using the function `red` applied to the color. Since internally the actual value is stored in a special format, we choose to convert it to a standard floating-point number using the `Float64` function:
Float64(Images.red(x))
#'

import Statistics
STAT = Statistics

funcs = [Images.red, Images.green, Images.blue]
imgs = [apple, banana]
A = [ STAT.mean(float.( f.(img) )) for f = funcs, img = imgs ]
#'

import Plots
import Plotly
#Plots.gr()
Plots.plotly()
v = float.(Images.green.(apple[:]))
pt = Plots.histogram(v, color="red", label="apple", normalize=true, nbins=25)
v = float.(Images.green.(banana[:]))
Plots.histogram!(pt, v, color="yellow", label="banana", normalize=true, nbins=25)
#'
pixel = apple[40,60]

red = Float64(Images.red(pixel))
green = Float64(Images.green(pixel))
blue = Float64(Images.blue(pixel))
println("The RGB values are ($red, $green, $blue)")

pixel = apple[1,1]

red = Float64(Images.red(pixel))
green = Float64(Images.green(pixel))
blue = Float64(Images.blue(pixel))
println("The RGB values are ($red, $green, $blue)")

R_of_apple = Float64.(Images.red.(apple))
STAT.mean(R_of_apple)
#'

R_of_banana = Float64.(Images.red.(banana))
STAT.mean(R_of_banana)
#'

G_of_banana = Float64.(Images.green.(banana))
STAT.mean(G_of_banana)

Ans = Plots.histogram(R_of_apple[:],color=:red,label="redness in the apple")

#' Julia's [mathematical standard library](https://docs.julialang.org/en/stable/stdlib/math/#Mathematics-1) has many mathematical functions built in. One of them is the `mean` function, which computes the average value. If we apply this to our apple:

STAT.mean(R_of_apple)

#' ## A quick riddle
#' Here's a quick riddle. Let's check the average value of red in the image of the banana.
R_of_banana = Float64.(Images.red.(banana))
STAT.mean(R_of_banana)
#' Oh no, that's more red than our apple? This isn't a mistake and is actually true! Before you move onto the next exercise, examine the images of the apple and the banana very carefully and see if you can explain why this is expected.

#' #### Exercise 1
#'
#' What is the average value of blue in the banana?
B_of_banana = Float64.(Images.blue.(banana))
STAT.mean(B_of_banana)

#' #### Exercise 2
#'
#' Does the banana have more blue or more green?
G_of_banana = Float64.(Images.green.(banana))
STAT.mean(G_of_banana)
#' which gives approximately `0.88`. The banana has more green on average.


text:
using LazyArtifacts
using DataFrames
using SQLite

db = SQLite.DB(joinpath(artifact"synpuf-10p", "synpuf-10p.sqlite"));
res = DBInterface.execute(db, "SELECT * FROM person");
println(DataFrame(res))


text:
# Function for reading in numerical parameters for basin simulations
function read_params(f_name)
    f = open(f_name, "r")
    params = []
    while ! eof(f)
        s = readline(f)
        if s[1] != '#'
            push!(params, split(s, '=')[2])
            flush(stdout)
        end
    end
    close(f)
    p = parse(Int64, params[1])
    T = parse(Float64, params[2])
    N = parse(Int64, params[3])
    Lw = parse(Float64, params[4])
    r̂ = parse(Float64, params[5])
    l = parse(Float64, params[6])
    b_depth = parse(Float64, params[7])
    dynamic_flag = parse(Int64,params[8])
    d_to_s = parse(Float64, params[9])
    dt_scale = parse(Float64, params[10])
    ic_file = params[11]
    ic_t_file = params[12]
    Dc = parse(Float64, params[13])
    B_on = parse(Int64, params[14])
    dir_out = params[15]
    volume_plots = parse(Int64, params[16])
    cycle_flag = parse(Int64, params[17])
    num_cycles = parse(Int64, params[18])
    intime_plotting = parse(Int64, params[19])
    μ_in = parse(Float64, params[20])
    return p, T, N, Lw, r̂, l, b_depth, dynamic_flag, d_to_s, dt_scale, ic_file, ic_t_file, Dc, B_on, dir_out, volume_plots, cycle_flag, num_cycles, intime_plotting, μ_in
end


text:
# Find the largest possible axes-aligned rectangle of true values within an array.
function largest_rect(X::AbstractArray{Bool})
	m,n = size(X)
	H = count.(X[1,:])
	col,height,areamax = largest_rect(H)
	rowrange = 1:1
	colrange = col .+ (0:height-1)
	for i in 2:m
		# in each column, heights of consecutive ones from row 1 to i
		H = [ X[i,j]*(1+H[j]) for j in 1:n ]
		newcol,newheight,newarea = largest_rect(H)
		if newarea > areamax
			areamax = newarea
			rowrange = i-newheight+1:i
			colrange = newcol .+ (0:newarea÷newheight-1)
		end
	end
	return rowrange,colrange,areamax

end

# Find the largest rectangle underneath a histogram-count vector.
function largest_rect(x::AbstractVector{T} where T <: Number)
    maxarea = 0
    x = [-1;x;-1]
    stack = [1]  # index into x

	idx = 0
	height = 0
    for i in 1:length(x)
        while x[i] < x[stack[end]]
			h = x[pop!(stack)]
			area = h*(i-stack[end]-1)
			if area > maxarea
				idx,height = stack[end],h
				maxarea = area
			end
		end
		push!(stack,i)
	end
	return idx,height,maxarea
end

# Grow a rectangular region to include all the adjacent pixels with similar value.
function grow_rectangle(X::AbstractMatrix,irange,jrange)
	m,n = size(X)
	rect = Array(CartesianIndices((irange,jrange)))
	isempty(rect) && (return rect)

	# Go far outside to get an idea of the intensity of non-belonging pixels in the vicinity.
	h,w = length(irange),length(jrange)
	mid = round.(Int,(median(irange),median(jrange)))
	ii = clamp.(mid[1]-h:mid[1]+h,1,m)
	jj = clamp.(mid[2]-w:mid[2]+w,1,n)
	outerval = median( [ vec(X[ii[[1,end]],jj]); vec(X[ii,jj[[1,end]]]) ] )

	# Define region and average value in it.
	region = falses(size(X))
	region[irange,jrange] .= true
	checked = copy(region)

	# Include those more like the average of the interior than the sampled exterior.
	avg = mean(X[region])
	thresh = max(0.2,0.5*(avg-outerval))
	nbrs = CartesianIndex.([(-1,0),(1,0),(0,-1),(0,1)])

	for kk = 1:1000
		added = 0
		for nbr in nbrs, r in findall(region)
			newpt = r + nbr
			if checkbounds(Bool,X,newpt) && !checked[newpt]
				if abs(X[newpt] - avg) < thresh
					region[newpt] = true
					added += 1
				end
				checked[newpt] = true
			end
		end
		#@show added
		(added == 0) && break
	end
	return findall(region)
end

"""
	findpurkinje(img,thresh=0.5)
Return indices of pixels believed to be in the purkinje of the image `img`. Decrease `thresh` toward zero to include more pixels, or increase to one to include fewer. Returns a vector of `CartesianIndex`.
"""
function findpurkinje(img::AbstractMatrix{T} where T<:Colorant,thresh=0.5;options=get_defaults(),rectangle=false)
	X = options.purkinje_channel.(img)
	idx = findpurkinje(X,thresh;options,rectangle=rectangle)
	return idx
end

function findpurkinje(X::AbstractMatrix{T} where T<:Number,thresh=0.5;options=get_defaults(),rectangle=false)
	m,n = size(X)
	θ,area = 0.95,0
	iran = jran = cen = NaN
	B = X.>θ
	ctr = 0
	removed = CartesianIndex{2}[]
	while θ >= thresh
		iran,jran,area = largest_rect(B)
		@debug "iran = $iran, jran = $jran"
		# If area is too small, there aren't enough pixels at this brightness.
		if (ctr > 4) || (area ≤ m*n*options.purkinje_minarea)
			θ -= 0.05
			ctr = 0
			B = X.>θ
			#B[removed] .= false
			@debug "θ = $θ"
		else
			h,w = (iran[end]-iran[1]+1,jran[end]-jran[1]+1)
			# Must have appropriate aspect ratio. Bright lines at eyelids can give short, skinny rectangles.
			if h > 0.8w
				# Must have contrast with neighboring pixels.
				cen = round.(Int,(median(iran),median(jran)))
				# Expand to larger box
				rows = clamp.(cen[1]-h:cen[1]+h,1,m)
				cols = clamp.(cen[2]-w:cen[2]+w,1,n)
				inner = sum(X[iran,jran])
				outer = sum(X[rows,cols]) - inner
				inneravg = inner/area
				outeravg = outer/((2h+1)*(2w+1)-area)
				@debug "inner = $inneravg, outer = $outeravg"
				if inneravg - outeravg > 0.25
					# we've got it
					break
				end
			end
			# Remove this region from consideration.
			bad = CartesianIndices((iran,jran))
			#bad = grow_rectangle(B,iran,jran)
			B[bad] .= false
			append!(removed,bad)
			@debug length(removed)
			ctr += 1
		end
	end
	return rectangle ? CartesianIndices((iran,jran)) : grow_rectangle(X,iran,jran)
end

"""
	findpeaks(index,sz,dim)
Look for peaks in the distribution of a vector `index` of indexes. Tuple `sz` is the size of the image that the indexes are taken from, and `dim` is the dimension (1=rows,2=columns) to operate in. Returns a pair of real values at the peaks of the distributions in the two halves of the range of the index.
"""
function findpeaks(idx,sz,dim)
	# look for the peaks in a distribution of indices
	vals = getindex.(idx,dim)

	edges = LinRange(1,sz[dim]+1,101)
	hist = fit(Histogram,vals,edges)
	w = hist.weights

	## search left/right sides separately
	k = argmax(w[1:50])
	jleft = mean(edges[k:k+1])
	wleft = w[k]
	k = argmax(w[51:end])
	jright = mean(edges[50+k:50+k+1])
	wright = w[50+k]

	return (jleft,wleft),(jright,wright)
end

max_stars_repo_name: tobydriscoll/cornea-detection

max_issues_repo_name: tobydriscoll/cornea-detection

max_forks_repo_name: tobydriscoll/cornea-detection


text:
export
    solve!,
    CBSSolver,
    MetaAgentCBS_Solver

function solve! end
function low_level_search! end

"""
    `solve!(solver, args ...)`

    Run the algorithm represented by `solver` on an instance of a Multi-Agent
    Path-Finding problem.
"""
function solve!(solver, args...)
    throw(ArgumentError(string(
        "function CRCBS.solve!(solver::", typeof(solver),",...) not defined.",
        " You must explicitly override CRCBS.solve!() for solvers of type",
        typeof(solver)
        )))
end

export
    BiLevelPlanner,
    AbstractCBSSolver,
    low_level

abstract type BiLevelPlanner end
low_level(solver::BiLevelPlanner) = solver.low_level_planner
abstract type AbstractCBSSolver <: BiLevelPlanner end
function set_best_cost!(solver::AbstractCBSSolver,cost)
    set_best_cost!(get_logger(solver),cost)
    set_best_cost!(low_level(solver),cost)
end
function hard_reset_solver!(solver::AbstractCBSSolver)
    hard_reset_solver!(get_logger(solver))
    hard_reset_solver!(low_level(solver))
end


text:
@testset "aggregation" begin
    v = rand(5)
    aggregate(v, mav) ≈ mean(v)
    aggregate(v, TruePositive()) ≈ sum(v)
    aggregate(v, rms) ≈ sqrt(mean(v.^2))
end

@testset "metadata" begin
    measures()
    measures(m -> m.target_scitype <: AbstractVector{<:Finite} &&
                  m.supports_weights)
end
true


text:
# ---
# title: 998. Maximum Binary Tree II
# id: problem998
# author: Tian Jun
# date: 2020-10-31
# difficulty: Medium
# categories: Tree
# link: <https://leetcode.com/problems/maximum-binary-tree-ii/description/>
# hidden: true
# ---
#
# We are given the `root` node of a _maximum tree:_ a tree where every node has
# a value greater than any other value in its subtree.
#
# Just as in the [previous problem](https://leetcode.comhttps://leetcode.com/problems/maximum-
# binary-tree/), the given tree was constructed from an list `A` (`root =
# Construct(A)`) recursively with the following `Construct(A)` routine:
#
#   * If `A` is empty, return `null`.
#   * Otherwise, let `A[i]` be the largest element of `A`.  Create a `root` node with value `A[i]`.
#   * The left child of `root` will be `Construct([A[0], A[1], ..., A[i-1]])`
#   * The right child of `root` will be `Construct([A[i+1], A[i+2], ..., A[A.length - 1]])`
#   * Return `root`.
#
# Note that we were not given A directly, only a root node `root =
# Construct(A)`.
#
# Suppose `B` is a copy of `A` with the value `val` appended to it.  It is
# guaranteed that `B` has unique values.
#
# Return `Construct(B)`.
#
#
#
# **Example 1:**
#
# **![](https://assets.leetcode.com/uploads/2019/02/21/maximum-binary-
# tree-1-1.png)![](https://assets.leetcode.com/uploads/2019/02/21/maximum-
# binary-tree-1-2.png)**
#
#
#
#     Input: root = [4,1,3,null,null,2], val = 5
#     Output: [5,4,null,1,3,null,null,2]
#     Explanation: A = [1,4,2,3], B = [1,4,2,3,5]
#
#
# **Example 2:**
#
# **![](https://assets.leetcode.com/uploads/2019/02/21/maximum-binary-
# tree-2-1.png)![](https://assets.leetcode.com/uploads/2019/02/21/maximum-
# binary-tree-2-2.png)**
#
#
#
#     Input: root = [5,2,4,null,1], val = 3
#     Output: [5,2,4,null,1,null,3]
#     Explanation: A = [2,1,5,4], B = [2,1,5,4,3]
#
#
# **Example 3:**
#
# **![](https://assets.leetcode.com/uploads/2019/02/21/maximum-binary-
# tree-3-1.png)![](https://assets.leetcode.com/uploads/2019/02/21/maximum-
# binary-tree-3-2.png)**
#
#
#
#     Input: root = [5,2,3,null,1], val = 4
#     Output: [5,2,4,null,1,3]
#     Explanation: A = [2,1,5,3], B = [2,1,5,3,4]
#
#
#
#
# **Constraints:**
#
#   * `1 <= B.length <= 100`
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end


text:
function wall_equation(cur_reactor::AbstractReactor)
  cur_G_T = K_WL(cur_reactor)

  cur_G_T *= cur_reactor.sigma_v ^ (1/3)

  cur_equation = Equation(1, 0, -2/3, cur_G_T)

  cur_equation
end


text:
using Isaac
using Test
using LinearAlgebra, SparseArrays, Random

print_tf(::Test.Pass) = printstyled("+", bold=true, color=:green)
print_tf(::Test.Fail) = printstyled("-", bold=true, color=:red)
print_tf(::Tuple{Test.Error,Bool}) = printstyled("x", bold=true, color=:magenta)

Random.seed!(1)

isCI = haskey(ENV, "CI")
notCI = !isCI

println("Running tests for `Isaac.jl`")
@testset "Isaac" begin

include("testproblems.jl")
include("test_nsoli.jl")
include("test_dlanczos.jl")
# include("test_saddle1.jl")
# include("test_minim.jl")

# if notCI
#     include("performance.jl")
# end

end