text:
# startup2.jl
# Merging mesh movement and FreeFormDeformation

push!(LOAD_PATH, "../src/")
push!(LOAD_PATH, joinpath(Pkg.dir("PumiInterface"), "src"))
push!(LOAD_PATH, joinpath(Pkg.dir("SummationByParts"), "src"))
push!(LOAD_PATH, joinpath(Pkg.dir("MeshMovement"), "src"))
push!(LOAD_PATH, joinpath(Pkg.dir("PDESolver"), "src/Utils"))

gc()

using PdePumiInterface
using SummationByParts
using ODLCommonTools
using ArrayViews
using Utils
using MPI
using MeshMovement
using FreeFormDeformation

# MPI Declarations
comm = MPI.COMM_WORLD
comm_world = MPI.MPI_COMM_WORLD
comm_self = MPI.COMM_SELF
my_rank = MPI.Comm_rank(comm)
comm_size = MPI.Comm_size(comm)

opts = PdePumiInterface.get_defaults()
# 2D mesh
opts["order"] = 1
opts["dimensions"] = 2
opts["use_DG"] = true
opts["operator_type"] = "SBPOmega"
opts["dmg_name"] = "./mesh_files/2DAirfoil_chord1.dmg"
opts["smb_name"] = "./mesh_files/2DAirfoil_chord1.smb"
opts["numBC"] = 2

# For boxInBox
#=
opts["BC1"] = [17,20,23,26]
opts["BC1_name"] = "FarField"
opts["BC2"] = [5,8,11,14]
opts["BC2_name"] = "Airfoil"
=#
# For 2DAirfoil
opts["BC1"] = [8]
opts["BC1_name"] = "FarField"
opts["BC2"] = [5]
opts["BC2_name"] = "Airfoil"

opts["coloring_distance"] = 2 # 0 For CG Mesh 2 for DG Mesh
opts["jac_type"] = 2

sbp, mesh, pmesh, Tsol, Tres, Tmsh, mesh_time = createMeshAndOperator(opts, 1)

geom_faces = opts["BC2"]

# Free Form deformation parameters
ndim = 2
order = [4,4,2]  # Order of B-splines in the 3 directions
nControlPts = [6,6,2]

ffd_map = PumiMapping{Tmsh}(ndim, order, nControlPts, mesh, full_geom=false,
geom_faces=geom_faces)

calcKnot(ffd_map)

# Create Bounding box
offset = [0., 0., 0.5] # No offset in the X & Y direction
ffd_box = PumiBoundingBox{Tmsh}(ffd_map, mesh, sbp, offset)

# Control points
controlPoint(ffd_map, ffd_box)
#= # test basis function evaluation
U = [0.0,0.0,0.0,0.0,0.3333333333333333,0.6666666666666666,1.0,1.0,1.0,1.0]
u = 0.5
span = FreeFormDeformation.findSpan(u,U,order[1],nControlPts[1])
println("span = $span")
N = zeros(order[1])
FreeFormDeformation.basisFunctions(U, order[1], u, span, N)
println("N =\n", N)
=#

# Populate map.xi
calcParametricMappingNonlinear(ffd_map, ffd_box, mesh, geom_faces)

# Prep MeshWarping
volNodes = zeros(Tmsh, 3, mesh.numVert)
MPI.Barrier(comm)
println("size mesh.vert_coords = $(size(mesh.vert_coords))")
for i = 1:mesh.numEl
  for j = 1:size(mesh.vert_coords,2)
    # Get the vertex numbering on the portion of mesh owned by the processor
    local_vertnum = mesh.element_vertnums[j,i]
    volNodes[1:2,local_vertnum] = mesh.vert_coords[:,j,i]# mesh.element_vertnums
    #println("Rank = $my_rank, volNodes[:, $local_vertnum] = $(volNodes[:,local_vertnum])")
  end
end
MPI.Barrier(comm)

# Get the bool array of all the surface vertex coordinates
surfaceVtx = zeros(Int32, mesh.numVert)
for (bindex, bndry) in enumerate(mesh.bndryfaces)
  vtx_arr = mesh.topo.face_verts[:, bndry.face]
  for j = 1:length(vtx_arr)
    local_vertnum = mesh.element_vertnums[vtx_arr[j], bndry.element]
    surfaceVtx[local_vertnum] = 1
  end
end

# Prepare the wall coordinates array for mesh warping
nwall_faces = zeros(Int,length(geom_faces))
vtx_per_face = mesh.dim # only true for simplex elements
for itr = 1:length(geom_faces)
  geom_face_number = geom_faces[itr]
  # get the boundary array associated with the geometric edge
  itr2 = 0
  for itr2 = 1:mesh.numBC
    if findfirst(mesh.bndry_geo_nums[itr2],geom_face_number) > 0
      break
    end
  end
  start_index = mesh.bndry_offsets[itr2]
  end_index = mesh.bndry_offsets[itr2+1]
  idx_range = start_index:(end_index-1)
  bndry_facenums =sview(mesh.bndryfaces, idx_range) # faces on geometric edge i
  nwall_faces[itr] = length(bndry_facenums)
end      # End for itr = 1:length(geomfaces)

# Populate wallCoords
wallCoords = zeros(3, sum(nwall_faces)*vtx_per_face)
ctr = 1 # Counter for wall coordinates
for itr = 1:length(geom_faces)
  geom_face_number = geom_faces[itr]
  # get the boundary array associated with the geometric edge
  itr2 = 0
  for itr2 = 1:mesh.numBC
    if findfirst(mesh.bndry_geo_nums[itr2],geom_face_number) > 0
      break
    end
  end
  start_index = mesh.bndry_offsets[itr2]
  end_index = mesh.bndry_offsets[itr2+1]
  idx_range = start_index:(end_index-1)
  bndry_facenums =sview(mesh.bndryfaces, idx_range) # faces on geometric edge i
  nfaces = length(bndry_facenums)
  for i = 1:nfaces
    bndry_i = bndry_facenums[i]
    # get the local index of the vertices
    vtx_arr = mesh.topo.face_verts[:,bndry_i.face]
    for j = 1:length(vtx_arr)
      wallCoords[1:2,ctr] = mesh.vert_coords[:,vtx_arr[j],bndry_i.element]
      #println("wallCoords[:, $ctr] = $(wallCoords[:,ctr])")
      ctr += 1
    end  # End for j = 1:length(vtx_arr)
  end    # End for i = 1:nfaces
end

# Warp parameters
nFaceConnectivity = sum(nwall_faces)*vtx_per_face
# Get total number of vertices in the mesh
total_vtx = MPI.Allreduce(mesh.numVert, MPI.SUM, comm)

println("nFaceConnectivity = $nFaceConnectivity")
param = warpParam(nNodes=total_vtx, nLocalNodes=mesh.numVert, nFaceConnLocal=nFaceConnectivity,
                  nFaceSizesLocal=sum(nwall_faces))
# MPI Information for warping
mpiVar = MPIInfo(warp_comm_world=comm_world, myID=my_rank, nProc=comm_size)

# Populate entries of param
param.aExp = 2.
param.bExp = 2.
param.LdefFact = 1
param.alpha = 0.2
param.symmTol = 1e-4
param.errTol = 1e-4
param.cornerAngle = 30.
param.zeroCornerRotations = false
param.useRotations = false
param.evalMode = 5

# Symmetry Information
# A symmetry plane is defined a point and normal. It is necessary to ensure that
# Each eplane is independent of the other. Look at checkPlane by Gaetan in
# UnstructuredMesh.Py. The way this is done in the data structure is that, every
# plane is defined by a (3,2) array. [:,1] defines the point and [:,2] defines
# the unit normal vector.
symmetryPlanes = zeros(Float64,3,2)
faceConn = collect(Int32, 0:nFaceConnectivity-1)
flatWarpSurfPts = reshape(wallCoords, 3*size(wallCoords,2))
faceSizes = mesh.dim*ones(Int32,sum(nwall_faces))
# warp the mesh

# Initialize warping
initializeWarping(param, mpiVar, symmetryPlanes, volNodes, surfaceVtx,
                      flatWarpSurfPts, faceConn, faceSizes)
# New Surface Coordinates
#=
# Print control points
for k = 1:ffd_map.nctl[3]
  for j = 1:ffd_map.nctl[2]
    for i = 1:ffd_map.nctl[1]
      println("cp_xyz[:,$i,$j,$k] = ", ffd_map.cp_xyz[:,i,j,k])
    end
    println('\n')
  end
end

# Camber the shit
ffd_map.cp_xyz[2,1:2,1:4,:] -= 0.01*4
ffd_map.cp_xyz[2,7,1:4,:] -= 0.015*4
ffd_map.cp_xyz[2,4:6,5:7,:] += 0.01*4

# Translate control points along x & y by  5 units
ffd_map.cp_xyz[1,:,:,:] += 0.1*4
ffd_map.cp_xyz[2,:,:,:] += 0.5*4

# Rotation matrix
theta = -20*pi/180  # Rotate wall coordinates by 10 degrees
rotMat = [cos(theta) -sin(theta) 0
          sin(theta) cos(theta)  0
          0          0           1] # Rotation matrix

# Rotate the control points
for k = 1:ffd_map.nctl[3]
  for j = 1:ffd_map.nctl[2]
    for i = 1:ffd_map.nctl[1]
      ffd_map.cp_xyz[:,i,j,k] = rotMat*ffd_map.cp_xyz[:,i,j,k]
    end
  end
end
=#

evalSurface(ffd_map, mesh)

# Update wallCoords
ctr = 1 # Counter for wall coordinates
for itr = 1:length(geom_faces)
  geom_face_number = geom_faces[itr]
  # get the boundary array associated with the geometric edge
  itr2 = 0
  for itr2 = 1:mesh.numBC
    if findfirst(mesh.bndry_geo_nums[itr2],geom_face_number) > 0
      break
    end
  end
  start_index = mesh.bndry_offsets[itr2]
  end_index = mesh.bndry_offsets[itr2+1]
  idx_range = start_index:(end_index-1)
  bndry_facenums =sview(mesh.bndryfaces, idx_range) # faces on geometric edge i
  nfaces = length(bndry_facenums)
  for i = 1:nfaces
    bndry_i = bndry_facenums[i]
    # get the local index of the vertices
    vtx_arr = mesh.topo.face_verts[:,bndry_i.face]
    for j = 1:length(vtx_arr)
      wallCoords[1:2,ctr] = mesh.vert_coords[:,vtx_arr[j],bndry_i.element]
      #println("wallCoords[:, $ctr] = $(wallCoords[:,ctr])")
      ctr += 1
    end  # End for j = 1:length(vtx_arr)
  end    # End for i = 1:nfaces
end

# Warp The mesh
warpMesh(param, volNodes, flatWarpSurfPts)
MPI.Barrier(comm)

# Update all the mesh coordinates
println("type of volNodes = ", typeof(volNodes), " size = ", size(volNodes))
println("type of mesh.vert_coords = ", typeof(mesh.vert_coords), " size = ", size(mesh.vert_coords))
for i = 1:mesh.numEl
  for j = 1:size(mesh.vert_coords,2)
    # Get the vertex numbering on the portion of mesh owned by the processor
    local_vertnum = mesh.element_vertnums[j,i]
    mesh.vert_coords[:,j,i] = volNodes[1:2,local_vertnum] # mesh.element_vertnums
  end
end

for i = 1:mesh.numEl
  update_coords(mesh, i, mesh.vert_coords[:,:,i])
end
commit_coords(mesh, sbp)
writeVisFiles(mesh, "visRotation")


text:
########################## Sparse Bayesian Learning ############################
# sparse bayesian learning for linear regression
# IDEA: unification of all sbl algorithms via functional α_update
struct SparseBayesianLearning{T, AT, BT, GT} <: Update{T}
    AΣA::AT
    AΣb::BT
    Γ::GT
    # B # pre-allocate
    # Σ::ST
    function SparseBayesianLearning(A::AbstractMatrix, b::AbstractVector,
                        Σ::AbstractMatrix, Γ::Diagonal = (1.0I)(size(A, 2)))
        AΣ = A'inverse(Σ)
        AΣb = AΣ * b
        AΣA = AΣ * A
        T = eltype(A)
        new{T, typeof(AΣA), typeof(AΣb), typeof(Γ)}(AΣA, AΣb, Γ)
    end
end
const SBL = SparseBayesianLearning
# IDEA: prune indices? with SparseArrays
# for better conditioning
# C = (get_sqrt(Γ)*AΣA*get_sqrt(Γ) + I)
# B = inverse(get_sqrt(Γ)) * C * inverse(get_sqrt(Γ))
# CG could work well, because we have good initial x's between iterations
# for noiseless L0 norm minimization: B = get_sqrt(Γ)*pseudoinverse(A * get_sqrt(Γ))
function update!(S::SBL, x::AbstractVector = zeros(size(S.Γ, 2)))
    AΣA, AΣb, Γ = S.AΣA, S.AΣb, S.Γ
    B = AΣA + inv(Γ) # woodury makes sense if number of basis functions is larger than x
    B = factorize(B) # could use in-place cholesky!(B)
    B isa Number ? (x .= AΣb ./ B) : ldiv!(x, B, AΣb) # x = B \ AΣb, could generalize with functional solve! input
    # update rules
    # @. Γ.diag = x^2 + $diag($inverse(B)) + 1e-16 # provably convergent
    @. Γ.diag = x^2 / (1 - $diag($inverse(B)) / Γ.diag) + 1e-14 # good heuristic
    return x
end

sbl(A, b, σ²::Real) = sbl(A, b, σ²*I(length(b)))
function sbl(A::AbstractMatrix, b::AbstractVector, Σ::AbstractMatrix;
                        maxiter::Int = 128size(A, 2), min_change::Real = 1e-6)
    P = SBL(A, b, Σ)
    x = zeros(size(A, 2))
    γ_old = copy(P.Γ.diag)
    for i in 1:maxiter
        update!(P, x)
        if norm(γ_old - P.Γ.diag) < min_change
            break
        end
        γ_old .= P.Γ.diag
    end
    return x
end

# updating noise variance (if scalar)
# B⁻¹ = inverse(B)
# σ² = sum(abs2, b-A*x) / sum(i->B[i]/Γ[i], diagind(Γ)) # could replace with euclidean norm

##################### greedy marginal likelihood optimization ##################
# updates the prior variance which leads to the largest
# increase in the marginal likelihood
struct FastSparseBayesianLearning{T, AT<:AbstractMatrix{T}, BT, NT, A, CT} <: Update{T}
    A::AT
    b::BT
    Σ::NT
    α::A

    # temporaries
    S::A
    Q::A
    δ::A # holds difference in marginal likelihood for updating an atom
    C⁻¹::CT # TODO: use Woodbury type here
end
const FSBL = FastSparseBayesianLearning
function FSBL(A::AbstractMatrix, b::AbstractVector, Σ)
    α = fill(Inf, size(A, 2))
    δ = similar(α)
    Σ = factorize(Σ)
    ΣA = Σ \ A
    Q = ΣA' * b
    ΣA .*= A
    S = vec(sum(ΣA, dims = 1))
    C⁻¹ = Matrix(inv(Σ))
    FSBL(A, b, Σ, α, S, Q, δ, C⁻¹)
end

function FSBL(A::AbstractMatrix, b::AbstractVector, σ²::Real)
    FSBL(A, b, σ²*I(length(b)))
end

##################### Relevance Matching Pursuit (RMP_σ) #########################
# IDEA: could merge FSBL and RMPS types
struct RMPS{T, AT<:AbstractMatrix{T}, BT, NT, A, CT} <: Update{T}
    A::AT
    b::BT
    Σ::NT
    α::A

    # temporaries
    S::A
    Q::A
    δ::A # holds difference in marginal likelihood for updating an atom
    C⁻¹::CT # IDEA: use Woodbury type here
end
function RMPS(A::AbstractMatrix, b::AbstractVector, σ²::Real, α::AbstractVector = fill(Inf, size(A, 2)))
    RMPS(A, b, σ²*I(length(b)), α)
end

# α is the inverse of γ, the prior variances of x
function RMPS(A::AbstractMatrix, b::AbstractVector, Σ::AbstractMatOrFac,
                                    α::AbstractVector = fill(Inf, size(A, 2)))
    n, m = size(A)
    δ = similar(α)
    C = if all(isinf, α)
            Σ
        else
            i = .!isinf.(α)
            γ = inv.(α[i])
            Γ = Diagonal(γ)
            Ai = A[:, i]
            Woodbury(Σ, Ai, Γ, Ai')
        end
    C = factorize(C)
    CA = C \ A
    Q = CA' * b
    CA .*= A
    S = vec(sum(CA, dims = 1))
    C⁻¹ = Matrix(inv(C))
    RMPS(A, b, Σ, α, S, Q, δ, C⁻¹)
end

# TODO: this can be done more efficiently with P.C
function Base.getproperty(S::Union{FSBL, RMPS}, s::Symbol)
    if s == :x
        active = @. !isinf(S.α) # active basis patterns
        x = spzeros(eltype(S.A), size(S.A, 2))
        Ai = @view S.A[:, active]
        Σ⁻¹ = inverse(S.Σ)
        P = inverse(*(Ai', Σ⁻¹, Ai) + Diagonal(S.α[active]))
        x[active] .= P * (Ai' * (Σ⁻¹ * S.b))
        dropzeros!(x)
    else
        getfield(S, s)
    end
end

################################################################################
# optimizes prior variances, returns corresponding weight estimate
# min_increase is the minimum increase in the marginal likelihood below which
# the algorithm terminates
function fsbl(A, b, Σ; maxiter = 2size(A, 2), min_increase = 1e-6)
    P = FSBL(A, b, Σ)
    optimize!(P, maxiter = maxiter, min_increase = min_increase)
    P.x
end
# optimizes prior variances, returns FSBL object
function optimize!(P::FSBL; maxiter = 2size(A, 2), min_increase = 1e-2)
    for i in 1:maxiter
        update!(P, P.α)
        if maximum(P.δ) < min_increase
            break
        end
    end
    return P
end

function update!(P::FSBL, α::AbstractVector = P.α)
    i = sbl_index!(P) # index corresponding to highest increase in marginal likelihood
    active, relevant = isactive(P, i), isrelevant(P, i)
    if !active && relevant
        sbl_acquisition!(P, i)
    elseif active && !relevant
        sbl_deletion!(P, i)
    elseif active && relevant
        sbl_update!(P, i)
    end # else, no action required
    return α
end

# updates set of active basis functions
# greedy strategy: chooses update which increases marginal likelihood most (max δ)
function update_α!(α, δ, S, Q)
    k = argmax(δ)
    sk, qk = get_sq(S[k], Q[k], α[k])
    α[k] = optimal_α(sk, qk)
    return α
end

############################## basic helpers ###################################
isactive(P::FSBL, i::Int) = (P.α[i] < Inf)
function isrelevant(P::FSBL, i::Int)
    s, q = get_sq(P.S[i], P.Q[i], P.α[i])
    return s < q^2
end
# calculate small s, q from S, Q (see Tipping 2003)
# also known as sparsity and quality factors
function get_sq(S::Real, Q::Real, α::Real)
    (S, Q) .* (α < Inf ? α / (α-S) : 1)
end

function optimal_α(s::Real, q::Real)
    (s < q^2) ? s^2 / (q^2 - s) : Inf
end

function sbl_index!(P::FSBL)
    @. P.δ = delta(P.α, P.S, P.Q)
    i = argmax(P.δ)
end

# potential change in marginal likelihood, and best action
function delta(α::Real, S::Real, Q::Real)
   s, q = get_sq(S, Q, α)
   active = α < Inf
   relevant = s < q^2
   if !active && relevant # out of model and high quality factor
       return δ_add(S, Q)
   elseif active && !relevant # in model but not high enough quality
       return δ_delete(S, Q, α)
   elseif active && relevant # α re-estimation
       αn = optimal_α(s, q)
       return δ = δ_update(S, Q, α, αn)
   else # !active && !relevant
       return 0.
   end
end

############################# acquisition ######################################
function sbl_acquisition!(P::Union{FSBL, RMPS})
    @. P.δ = sbl_acquisition_value(P.α, P.S, P.Q)
    k = argmax(P.δ) # choosing maximum q/s value, corresponds to maximum increase in likelihood
    return sbl_acquisition!(P, k)
end
function sbl_acquisition!(P::Union{FSBL, RMPS}, i::Int)
    if P.δ[i] > 0 # only add if index i is beneficial
        si, qi = get_sq(P.S[i], P.Q[i], P.α[i]) # updating α
        P.α[i] = optimal_α(si, qi)
        γ = inv(P.α[i])
        update_SQC!(P, i, γ)
        return true
    end
    return false
end
function sbl_acquisition_value(α::Real, S::Real, Q::Real)
    s, q = get_sq(S, Q, α)
    active = α < Inf
    relevant = s < q^2
    return (!active && relevant) ? δ_add(S, Q) : 0.
end

function δ_add(S::Real, Q::Real)
   (Q^2 - S) / S + log(S) - log(Q^2)
end

function rmp_acquisition_value(α::Real, S::Real, Q::Real)
    s, q = get_sq(S, Q, α)
    active = α < Inf
    relevant = s < q^2
    if !isactive && isrelevant
        q^2 / s
    else
        0.
    end
end

################################ update ########################################
# delete index with maximum increase in likelihood
function sbl_update!(P::Union{FSBL, RMPS})
    @. P.δ = sbl_update_value(P.α, P.S, P.Q) # only necessary for active atoms
    i = argmax(P.δ)
    return sbl_update!(P, i)
end
# delete index i, if it increases the likelihood
function sbl_update!(P::Union{FSBL, RMPS}, i::Int)
    if P.δ[i] > 0
        si, qi = get_sq(P.S[i], P.Q[i], P.α[i])
        αs = optimal_α(si, qi)
        γ = inv(αs) - inv(P.α[i]) # difference between new γ and old γ
        update_SQC!(P, i, γ)
        P.α[i] = αs
        return P.δ[i]
    end
    return zero(eltype(P.δ[i]))
end

function sbl_update_value(α::Real, S::Real, Q::Real)
    α < Inf || return zero(α) # is active
    s, q = get_sq(S, Q, α)
    s < q^2 || return zero(α) # is relevant
    αn = optimal_α(s, q)
    return δ_update(S, Q, α, αn)
end
# potential change in marginal likelihood upon updating α to αn
function δ_update(S::Real, Q::Real, α::Real, αn::Real)
   d = 1/αn - 1/α
   Q^2 / (S + 1/d) - log(max(1 + S*d, 0))
end

################################ deletion ######################################
function sbl_deletion!(P::Union{FSBL, RMPS})
    @. P.δ = sbl_deletion_value(P.α, P.S, P.Q)
    i = argmax(P.δ)
    return sbl_deletion!(P, i)
end
function sbl_deletion!(P::Union{FSBL, RMPS}, i::Int)
    if P.δ[i] > 0
        si, qi = get_sq(P.S[i], P.Q[i], P.α[i]) # updating α
        γ = -inv(P.α[i]) # difference between new γ (0) and old γ (1/α[k])
        P.α[i] = optimal_α(si, qi) # setting to Inf
        update_SQC!(P, i, γ)
        return true
    end
    return false
end
function sbl_deletion_value(α::Real, S::Real, Q::Real)
    α < Inf || return zero(α) # is active
    s, q = get_sq(S, Q, α)
    s > q^2 || return zero(α) # is irrelevant
    return δ_delete(S, Q, α)
end
# potential change in marginal likelihood upon setting α to Inf
function δ_delete(S::Real, Q::Real, α::Real)
   Q^2 / (S - α) - log(1 - (S/α)) # delete k
end

######################## computing and updating S, Q, C⁻¹ ######################
get_SQ!(P::RMPS) = get_SQ!(P.S, P.Q, P.A, P.b, P.C⁻¹)
function get_SQ!(S, Q, A, b, C⁻¹) # try to avoid
    for k in 1:size(A, 2)
        Ak = @view A[:, k]
        S[k] = dot(Ak, C⁻¹, Ak)
        Q[k] = dot(Ak, C⁻¹, b)
    end
    S, Q
end

# get kernel matrix C = Σ + A * Γ * A'
function get_C_inverse(isactive, A, α, Σ)
    if any(isactive)
        Ai = @view A[:, isactive]
        Γ = Diagonal(inv.(α[isactive]))
        C = Woodbury(Σ, Ai, Γ, Ai')
        C⁻¹ = inverse(factorize(C))
    else
        C⁻¹ = inverse(Σ)
    end
end

# update after rank one modification
function update_C_inverse!(P::Union{FSBL, RMPS}, i::Int, α::Real)
    a = @view P.A[:, i]
    v = P.C⁻¹ * a
    @. P.C⁻¹ -= v * v' / (α + P.S[i]) # update
end

# updates C⁻¹, S, Q in P corresponding to rank one correction of C
# with ith atom with coefficient γ: C + γ * ai * ai'
function update_SQC!(P::Union{FSBL, RMPS}, i::Int, γ::Real)
    α = inv(γ)
    a = P.A[:, i]
    v = P.C⁻¹ * a
    update_C_inverse!(P, i, α)
    Av = P.A' * v
    si, qi = P.S[i], P.Q[i]
    @. P.S -= Av^2 / (α + si)
    @. P.Q -= Av * qi / (α + si)
    return P.S, P.Q
end

################################################################################
# Σ is the noise variance (scalar or matrix)
# maxiter is the maximum number of outer iterations
# maxiter_acquisition is the maximum number of iterations per aquisition stage
# maxiter_deletion is the maximum number of iterations per deletion stage
# min_increase is the minimum increase in marginal likelihood, below which the algorithm terminates
function rmps(A, b, Σ; maxiter::Int = size(A, 1), maxiter_acquisition::Int = size(A, 1),
                    maxiter_deletion::Int = size(A, 1), min_increase::Real = 1e-6) # in RMP experiments this was min_increase = 1e-2
    P = RMPS(A, b, Σ)
    optimize!(P, maxiter = maxiter, maxiter_acquisition = maxiter_acquisition,
                maxiter_deletion = maxiter_deletion, min_increase = min_increase)
    return P.x
end

function optimize!(P::RMPS; maxiter::Int = size(P.A, 1),
                            maxiter_acquisition::Int = size(P.A, 1),
                            maxiter_deletion::Int = size(P.A, 1),
                            min_increase::Real = 1e-6)
    A, b, Σ, S, Q = P.A, P.b, P.Σ, P.S, P.Q
    α = P.α
    α .= Inf
    old_α = copy(α)
    for i in 1:maxiter
        for _ in 1:maxiter_acquisition # acquisition stage
            rmp_acquisition!(P) || break # if we did not add any atom, we're done
        end
        old_α != α || break # while norm of α is still changing
        old_α .= α
        for j in 1:maxiter_deletion # update and deletion stage
            if rmp_deletion!(P)
                continue
            elseif rmp_update!(P) < min_increase
                break
            end
        end
        old_α != α || break # while norm of α is still changing
        old_α .= α
    end
    return P
end

# acquisition and update are identical to fast sbl
rmp_acquisition!(P::RMPS) = sbl_acquisition!(P)
rmp_update!(P::RMPS) = sbl_update!(P)

################################ deletion ######################################
function rmp_deletion!(P::RMPS)
    α, S, Q = P.α, P.S, P.Q
    @. P.δ = rmp_deletion_value(α, S, Q)
    min_δ, k = findmin(P.δ) # choosing minimum q/s value
    if min_δ < 1 # only delete if we are not doing damage
        sk, qk = get_sq(S[k], Q[k], α[k]) # updating α
        γ = -inv(α[k]) # difference between new γ (0) and old γ (1/α[k])
        α[k] = optimal_α(sk, qk) # setting to Inf
        update_SQC!(P, k, γ)
        return true
    end
    return false
end

# WARNING: lower is better here
function rmp_deletion_value(α::Real, S::Real, Q::Real)
    s, q = get_sq(S, Q, α)
    isactive = α < Inf
    isrelevant = s < q^2
    if isactive && !isrelevant # in model but not high enough quality
        δ = q^2 / s
    else
        Inf
    end
end

##################### optimization of noise variance ###########################
# optimizes noise variance in an outer loop above rmps
# a_σ², b_σ² are shape and scale parameter for Inverse Gamma prior on σ²
function rmps(A::AbstractMatrix, b::AbstractVector, optimize_σ²::Val{true},
                        σ²::Real = 1e-2, a_σ²::Real = 0, b_σ²::Real = 0;
                        maxiter::Int = 2size(A, 2), min_increase::Real = 1e-6,
                        maxouteriter::Int = 16, min_change::Real = 1e-12)
    α = fill(Inf, size(A, 2))
    for i in 1:maxouteriter
        P = RMPS(A, b, σ², α)
        optimize!(P, maxiter = maxiter, min_increase = min_increase)
        α = P.α
        σ²_new = estimate_σ²(P, a_σ², b_σ²)
        converged = abs(σ²_new - σ²) < min_change
        σ² = σ²_new
        if converged
            break
        end
    end
    x = RMPS(A, b, σ², α).x
    return x, σ²
end

function estimate_σ²(P::RMPS, a_σ²::Real = 0, b_σ²::Real = 0)
    estimate_σ²(P.A, P.x, P.b, inv.(P.α), a_σ², b_σ²)
end
# with inverse-gamma prior on σ² with shape and scale parameters a, b
function estimate_σ²(A::AbstractMatrix, x::AbstractVector, b::AbstractVector,
                    γ::AbstractVector, a_σ²::Real = 0, b_σ²::Real = 0)
    n = size(A, 1)
    (sum(abs2, b-A*x) + 2b_σ²) / (n - sum(γ) + 2a_σ²)
end

ext: jl


text:
#!/usr/bin/julia
#
# Functions for statistical analysis.
#
import Statistics#: std;

export get_stats, get_stats!, get_stats_str, get_stats_str!, get_csv_stats_str;

"""
    get_stats(quantity::Vector)

Compute basic statistics (mean, minimum and maximum) on the collection of data
in `quantity`.

See also: [`get_stats!`](@ref)
"""
function get_stats(quantity::Vector)
  av = sum(quantity)/length(quantity);
  max = maximum(quantity);
  min = minimum(quantity);
  stdv = Statistics.std(quantity, mean=av);

  return (av, max, min, stdv);
end
"""
    get_stats!(quantity::Vector)

Remove zero elements from `quantity` then return `get_stats(quantity)`.

See also: [`get_stats`](@ref)
"""
function get_stats!(quantity::Vector)
  filter!(v->(v!=0), quantity);
  return get_stats(quantity);
end

"""
    get_stats_str(quantity::Vector, name::String="quantity"; long::Bool=true)

Construct a string containing a human readable presentation of the basic
statistics returned from `get_stats(quantity)` using the label `name`.

Passing `long=false` prints the statistics on a single (line feed terminated)
line.

See also: [`get_csv_stats_str`](@ref), [`get_stats_str!`](@ref), [`get_stats`](@ref)
"""
function get_stats_str(quantity::Vector, name::String="quantity"; long::Bool=true)
  return get_stats_str(get_stats(quantity)..., name, long=long);
end
function get_stats_str(av::Number, max::Number, min::Number, stdv::Number,
                       name::String="quantity"; long::Bool=true)
  # Construct string.
  if long
    return string("Average ", name, " = ", av, '\n',
                  "Standard deviation of ", name, " = ", stdv, '\n',
                  "Maximum ", name, " = ", max, '\n',
                  "Minimum ", name, " = ", min, '\n');
  else
    return string(uppercasefirst(name), " statistics:\n",
                  "μ=", av, ", σ=", stdv, ", max=", max, ", min=", min, '\n');
  end
end

"""
    get_stats_str!(quantity::Vector, name::String="quantity"; long::Bool=false)

Construct a string containing a human readable presentation of the basic
statistics returned from `get_stats!(quantity)` using the label `name`.

Passing `long=false` prints the statistics on a single (line feed terminated)
line.

See also: [`get_stats_str`](@ref), [`get_stats!`](@ref)
"""
function get_stats_str!(quantity::Vector, name::String="quantity";
                        long::Bool=true)
  return get_stats_str(get_stats!(quantity)..., name, long=long);
end

"""
    get_csv_stats_str(quantity::Vector, name::String="quantity";
                      label::Bool=true, newline::Bool=true)

Construct a string containing CSV values for the basic statistics returned from
`get_stats(quantity)` using the label `name` as the first CSV column entry.

See also: [`get_stats_str`](@ref), [`get_stats_str!`](@ref), [`get_stats`](@ref)
"""
get_csv_stats_str(quantity::Vector, name::String="quantity";
                  label::Bool=true, newline::Bool=true) =
  get_csv_stats_str(get_stats(quantity)..., name, label=label, newline=newline);
function get_csv_stats_str(av::Number, max::Number, min::Number, stdv::Number,
                           name::String="quantity"; label::Bool=true, newline::Bool=true)
  # Construct string.
  return string((label ? string(name,",") : ""),
                av,",",stdv,",",max,",",min,(newline ? '\n' : ','));
end
function get_csv_stats_str(av::T, max::T, min::T, stdv::T, name::String="quantity";
                           label::Bool=true, newline::Bool=true) where
  {T<:Unitful.Quantity}
  # Construct string.
  return string((label ? string(name," / ",unit(T),",") : ""),
                ustrip(av),",",ustrip(stdv),",", ustrip(max),",",ustrip(min),
                (newline ? '\n' : ','));
end


text:
using Fairness
using MLJBase, MLJModels
using Test

@testset "fair_tensor" begin
    include("fair_tensor.jl")
end

@testset "measures" begin
    include("measures/measures.jl")
end

@testset "dataset_macros" begin
    include("datasets/datasets.jl")
    include("datasets/synthetic.jl")
end

@testset "algorithms" begin
    include("algorithms/algorithms.jl")
end

@testset "benchmarking" begin
    include("benchmark/benchmark.jl")
end


text:
using IntervalWavelets
using Base.Test

include("Quad.jl")
include("DyadicPoints.jl")
include("Daubechies.jl")
include("BoundaryDaubechies.jl")
include("Reconstruction.jl")


text:
using MetabolicFluxAnalysis

# implements the volume computation algorithm presented in
# I. Z. Emiris and V. Fisikopoulos , 2014 , Efﬁcient random-walk methods for approximating polytope volume.
# currently, it uses slower (uniform direction) hr sampling instead of the
# proposed chr sampling. This will be implemented soon.
function vol_mmc_reversed( G , h , num_samples , n_steps_initial , n_steps_rep , rho )

    d = size(G,2)

    # num iterations to do, such that rho^m >= d^d
    m = Integer( ceil( ( d * log(d) ) / log(rho) ) )

    # construct the radii that we have to do in the rounded polytope..
    # i.e. all from 1.0 to
    r_m = rho.^( (0:m) / d )

    print("r_m: $(r_m)")

    # round polytope:
    if(true)
        #tic()
        # compute MVE
        # E = {L^(-1)*x + f: norm(x) <= 1}
        (f,L) = mve_run(G, h)
        # rescale to isotropic ball
        h = h - G*f
        G = (G/L)/sqrt(d+2)
        #time_rounding = toc()
    end

    # Generate initial samples:
    #tic()
    # warm-up samples..
    print("\nWarmup.. \n")
    #switch(hr_conf.method)
    #    case {'HitAndRun'}
    X0 = sample_unif_hr( MetabolicFluxAnalysis.Polyhedron(G,h[:,:]) , zeros(d,num_samples) ,  n_steps_initial )
    #    case 'CoordHR',
            #%X0 = fullCoordHR(G,h,100*r_m(end),zeros(d,1),num_samples,n_steps_rep);
            #X0 = fullCoordHR(G,h,100*r_m(end),zeros(d,1),num_samples,n_steps_initial_decorr);
    #end
    #time_initialsamples = toc()

    ## Run over the reverse mmc iterations..
    X = X0

    factors = ones(Integer(m)) #Vector{Float64}(m)
    n_out_stored = zeros(Integer(m))#Vector{Int64}(m)
    n_resampled_tot = 0

    #tic()
    # initial radius:
    for zi=m:-1:2
        # ball radius of iteration:
        r_old = sqrt(d+2) * r_m[zi]
        r_new = sqrt(d+2) * r_m[zi-1]

        idx_out = sqrt.( sum( X.^2 , 1 ) ) .> r_new
        n_out   = sum(idx_out)
        ratio   = num_samples ./ (num_samples-sum(idx_out))

        n_out_stored[zi] = n_out
        # store result
        factors[zi] = (num_samples/(num_samples-n_out))
        n_resampled_tot = n_resampled_tot + n_out
        # Remove samples that are out:
        X(:,idx_out) = []

        # Resample:
        if( (n_out>0) && (zi > 2) )
            # select random samples which will be decorrelated..
            idx_rep = randsample( size(X,2) , n_out ) # should work here
            X_rep = X[ : , idx_rep ]

            # Run Hit-and-run steps:
            #switch(hr_conf.method),
            #    case {'HitAndRun','CoordinateHR'},
            #X_rep = XUniformHitAndRun_SingleCore( X_rep , struct('A',G,'b',h) , struct('L',eye(d),'f',zeros(d,1),'r',r_new,'L_is_eye',1) , n_steps_rep , hr_conf , hr_output_conf  );
            #X_rep = hrstep!( G , h , r_new , X_rep , n_steps_rep)
            X_rep = hrstep_cchr!!( G , h , r_new , X_rep , n_steps_rep)
            #    case 'MexCoordHR',
            #        X_rep = fullCoordHR(G,h,r_new,X_rep,1,n_steps_rep);
            #end
            #X[:,end+1:end+n_out] = X_rep;
            X = [ X  X_rep ]
        end
        # it/allit   factora factorb   tot_resampled
        #fprintf(' %04d / %04d   %d/%d    %d    \n',   m-zi+1,m  , num_samples,num_samples-n_out , n_resampled_tot );

        #if(total_debug), plotcolors = linspecer(m); plot(X(1,1:500),X(2,1:500),'.','Color',plotcolors(zi,:)); end
    end
    #time_mmc = toc()

    # original volume of ellipsoid corresopnding to the smallest ball:
    #lvol = (d/2)*log(pi) - gammaln(d/2+1) - sum(log(diag(L)))
    lvol = (d/2)*log(pi) - lgamma(d/2+1) - sum(log.(diag(L)))
    #%lvol = (d/2)*log(pi) - gammaln(d/2); %- sum(log(diag(L)));

    # multiply with the estimated factors..
    lvol = lvol + sum( log.(factors[2: length(factors) ] ) )

    logvol = lvol
    vol    = exp(lvol)
    #MoreInfo.time_Rounding       = time_rounding;
    #MoreInfo.time_InitialSamples = time_initialsamples;
    #MoreInfo.time_MMC            = time_mmc;
    #MoreInfo.time_Total          = toc(tic_start);

    #MoreInfo.factors             = factors(2:end);
    #MoreInfo.n_out               = n_out_stored(2:end);
    #MoreInfo.logvol_mve          = (d/2)*log(pi) - gammaln(d/2+1) - sum(log(diag(L)));
end


text:
# ============read and write segy file ====================
export convert,
       show,
       read_text_header,
       print_text_header,
       put_text_header,
       FileHeader,
       init_FileHeader,
       read_file_header,
       put_file_header,
       TraceHeader,
       init_TraceHeader,
       extract_traces_header,
       write_traces_header,
       read_traces_header,
       extract_field_trace_header,
       print_field_trace_header,
       read_segy_file,
       write_segy_file,
       write_one_trace,
       extract_sensors_location,
       create_acquisition_geometry,
       create_lookup_table,
       read_lookup_table,
       read_one_shot,
# ======Regular Sampled data format========================
       RegularSampleHeader,
       zeros,
       field_location,
       write_RSheader,
       write_RSdata,
       read_RSheader,
       read_RSdata,
       l2norm_rsf

include("segy.jl")
include("regular.jl")


text:
"""
    @fastfunc functionname[::String] polynomial[::Polynomial{true}]

    Defines a function with name `functionname` for `polynomial` for fast numerical evaluation.
    Note that the polynomial is transformed using `@fastmath`, which calls functions that
    may violate strict IEEE semantics.

    # Examples
    ```
    julia> using SphericalHarmonics

    julia> @polyvar x y z
    (x, y, z)

    julia> p = 15.0*x*y^2+7.5*x*z^13
    7.5xz¹³ + 15.0xy²

    julia> foo = @fastfunc p;

    julia> foo(1.0,2.0,3.0)
    1.19574825e7
    ```
    """
macro fastfunc(polynomial)
	return quote
		local polystr = string($(esc(polynomial)))
		local vars = string(tuple(variables($(esc(polynomial)))...))
		# insert "*" in between numbers and variables
		local regex = Regex("(\\d)" * replace(vars,", " => "|"))
		polystr = replace(polystr,regex => s"\1*\2")
		# insert "*" in between variables "xy" -> "x*y"
		local regex_xy = Regex(replace(vars,", " => "|")*replace(vars,", " => "|"))
		for i=1:length(variables($(esc(polynomial))))-1
			polystr = replace(polystr,regex_xy => s"\1*\2")
		end
		# replace "+ -" with "- "
		polystr = replace(polystr,r"([+])(\s)([-])" => s"\3 ")
		# create expression for function definition

		eval(Meta.parse(vars*" -> @fastmath "*polystr))
	end
end


text:
using Test
using ReadStatTables

using CategoricalArrays
using DataFrames
using Dates
using ReadStat
using ReadStatTables: _parse_usecols
using Tables

const tests = [
    "LabeledArrays",
    "table",
    "readstat"
]

printstyled("Running tests:\n", color=:blue, bold=true)

@time for test in tests
    include("$test.jl")
    println("\033[1m\033[32mPASSED\033[0m: $(test)")
end


text:
###########################################################################
#  Slice sampler (SliceSampler)
#
#  References:
#    - Neal R.M. Slice Sampling. The Annals of Statistics, 2003, 31 (3), pp 705–1031
#    - MacKay model.size. Information Theory, Inference, and Learning Algorithms. Cambridge University Press, 2003, section 29.7,
#      pp 374-378
#
#  A slice sampler is an adaptive step-size MCMC algorithm for continuous random variables, which only requires an
#  unnormalized density function as input.  It is a convenient alternative to Metropolis because it often gives good
#  results with very little tuning, and works in the case that different parts of the distribution require different
#  step sizes.  However, if the distribution's widths don't vary too much, and there's a good initialization and
#  proposal, Metropolis may be more efficient. This slice sampler is univariate; for a multivariate problem it just
#  updates variables one at a time.  Thus it suffers when variables are correlated.
#
#  Parameters :
#    - widths: Step sizes for initially expanding the slice (model.size-dim vector).
#    - stepout: Protects against the case of passing in small widths.
#
#  This is a port of Iain Murray's implementation, who kindly accepted to share his code at
#  http://homepages.inf.ed.ac.uk/imurray2/teaching/09mlss/slice_sampler.m
#
###########################################################################

export SliceSampler

println("Loading SliceSampler(widths, stepout)")

###########################################################################
#                  SliceSampler type
###########################################################################

immutable SliceSampler <: MCMCSampler
  widths::Union(Vector{Float64}, Nothing)
  stepout::Bool

  function SliceSampler(widths::Union(Vector{Float64}, Nothing), stepout::Bool)
    # @assert x > 0 "widths should be > 0"
    new(widths, stepout)
  end
end
SliceSampler() = SliceSampler(nothing, true)
SliceSampler(widths::Union(Vector{Float64}, Nothing)) = SliceSampler(widths, true)
SliceSampler(stepout::Bool) = SliceSampler(nothing, stepout)
SliceSampler(; widths::Union(Vector{Float64}, Nothing)=nothing, stepout::Bool=true) = SliceSampler(widths, stepout)

###########################################################################
#                  SliceSampler task
###########################################################################

# SliceSampler sampling
function SamplerTask(model::MCMCModel, sampler::SliceSampler, runner::MCMCRunner)
  local state::Vector{Float64}, x_l::Vector{Float64}, x_r::Vector{Float64}, xprime::Vector{Float64}
  local logTarget::Float64, log_uprime::Float64, r::Float64
  local widths::Vector{Float64}
  if sampler.widths == nothing
    widths = ones(model.size)
  else
    @assert length(sampler.widths)==model.size "Length of step sizes in widths must be equal to model size"
    widths = sampler.widths
  end

  state = copy(model.init)
  logTarget = model.eval(state)

  i = 1
  while true
    for dd = 1:model.size
      log_uprime = log(rand()) + logTarget
      x_l  = copy(state)
      x_r  = copy(state)
      xprime = copy(state)

      # Create a horizontal interval (x_l, x_r) enclosing xx
      r = rand()
      x_l[dd] = state[dd] - r*widths[dd]
      x_r[dd] = state[dd] + (1-r)*widths[dd]
      if sampler.stepout
        while model.eval(x_l) > log_uprime
          x_l[dd] -= widths[dd]
        end
        while model.eval(x_r) > log_uprime
          x_r[dd] += widths[dd]
        end
      end

      # Inner loop: propose xprimes and shrink interval until good one is found
      while true
        xprime[dd] = rand() * (x_r[dd] - x_l[dd]) + x_l[dd];
        logTarget = model.eval(xprime)
        if logTarget > log_uprime
          break
        else
          if (xprime[dd] > state[dd])
            x_r[dd] = xprime[dd];
          elseif (xprime[dd] < state[dd])
            x_l[dd] = xprime[dd];
          else
            @assert false "BUG, shrunk to current position and still not acceptable";
          end
        end
      end

      state[dd] = xprime[dd]
    end

    ms = MCMCSample(state, logTarget, fill(NaN, model.size), NaN)
    produce(ms)

    i += 1
  end
end


text:
using Surrogates

#1D
x = [1.0,2.0,3.0]
y = [3.0,5.0,7.0]
lb = 0.0
ub = 5.0
my_wend = Wendland(x,y,lb,ub)
add_point!(my_wend,0.5,4.0)
val = my_wend(0.5)

#ND
lb = [0.0,0.0]
ub = [4.0,4.0]
x = sample(5,lb,ub,SobolSample())
f = x -> x[1]+x[2]
y = f.(x)
my_wend_ND = Wendland(x,y,lb,ub)
est = my_wend_ND((1.0,2.0))
add_point!(my_wend_ND,(3.0,3.5),4.0)
add_point!(my_wend_ND,[(9.0,10.0),(12.0,13.0)],[10.0,11.0])
#todo


text:
mutable struct MixedGrid{dim,C,T<:Real} <: JuAFEM.AbstractGrid
    cells::Vector{C}
    nodes::Vector{Node{dim,T}}
    # Sets
    cellsets::Dict{String,Set{Int}}
    nodesets::Dict{String,Set{Int}}
    facesets::Dict{String,Set{Tuple{Int,Int}}} # TODO: This could be Set{FaceIndex} which could result in nicer use later
    # Boundary matrix (faces per cell × cell)
    boundary_matrix::SparseMatrixCSC{Bool,Int}
end

function MixedGrid(
    cells::Vector{C},
    nodes::Vector{Node{dim,T}};
    cellsets::Dict{String,Set{Int}}=Dict{String,Set{Int}}(),
    nodesets::Dict{String,Set{Int}}=Dict{String,Set{Int}}(),
    facesets::Dict{String,Set{Tuple{Int,Int}}}=Dict{String,Set{Tuple{Int,Int}}}(),
          boundary_matrix::SparseMatrixCSC{Bool,Int}=spzeros(Bool, 0, 0)) where {dim,C,T}
    return MixedGrid(cells, nodes, cellsets, nodesets, facesets, boundary_matrix)
end


text:
using ClimateMachine

ClimateMachine.init()

using ClimateMachine.CartesianDomains

using CLIMAParameters: AbstractEarthParameterSet
struct DefaultParameters <: AbstractEarthParameterSet end

using ClimateMachine.Ocean:
    HydrostaticBoussinesqSuperModel, current_time, current_step, Δt

@testset "$(@__FILE__)" begin

    domain = RectangularDomain(
        Ne = (1, 1, 1),
        Np = 4,
        x = (-1, 1),
        y = (-1, 1),
        z = (-1, 0),
    )

    model = HydrostaticBoussinesqSuperModel(
        domain = domain,
        time_step = 0.1,
        parameters = DefaultParameters(),
    )

    @test model isa HydrostaticBoussinesqSuperModel
    @test Δt(model) == 0.1
    @test current_step(model) == 0
    @test current_time(model) == 0.0
end


text:
using Documenter, DiffEqFlux

makedocs(
    sitename = "DiffEqFlux.jl",
    authors="Chris Rackauckas et al.",
    clean = true,
    doctest = false,
    modules = [DiffEqFlux],

    format = Documenter.HTML(#analytics = "",
                             assets = ["assets/favicon.ico"],
                             canonical="https://diffeqflux.sciml.ai/stable/"),
    pages=[
        "Home" => "index.md",
        "Tutorials" => Any[
            "examples/LV-ODE.md",
            "examples/LV-stochastic.md",
            "examples/LV-Flux.md",
            "examples/LV-delay.md",
            "examples/NN-ODE.md",
            "examples/NeuralODE_Flux.md",
            "examples/Supervised-NN-ODE-MNIST.md",
            "examples/NN-CNF.md",
            "examples/local_minima.md",
            "examples/NN-SDE.md",
            "examples/NeuralOptimalControl.md",
            "examples/LV-Univ.md",
            "examples/SecondOrderNeural.md",
            "examples/LV-NN-Stiff.md",
            "examples/NewtonSecondOrderAdjoints.md",
            "examples/LV-Jump.md",
            "examples/universaldiffeq.md",
            "examples/minibatch.md"
        ],
        "Continuous Normalizing Flows Layer" => "CNFLayer.md",
        "Neural Differential Equation Layers" => "NeuralDELayers.md",
        "Use with Flux Chain and train!" => "Flux.md",
        "FastChain" => "FastChain.md",
        "GPUs" => "GPUs.md",
        "sciml_train" => "Scimltrain.md",
        "Benchmark" => "Benchmark.md"
    ]
)

deploydocs(
   repo = "github.com/SciML/DiffEqFlux.jl.git";
   push_preview = true
)


text:
gplfile3 = joinpath(dirname(@__FILE__),"fit3.gpl")

"""
    type InterpolationMethod <: AbstractRootsMethod
        dx::Float64
        maxiters::Int
        parallel::Bool
        verb::Int
        atol::Float64
        n::Int
    end

Type containg the parameters for root finding through interpolation.
`n`: number of points to interpolate at each iteration. `dx`: points interspacing
for the first interpolation.
If `dx<0` the value `abs(x0)/20` shall be used.

The default parameters are:

    InterpolationMathod(dx=-1, maxiters=20, parallel= false, verb=0, atol=1e-10, n=6)

"""
mutable struct InterpolationMethod <: AbstractRootsMethod
    dx::Float64
    maxiters::Int
    parallel::Bool
    verb::Int
    atol::Float64
    n::Int
end

InterpolationMethod(; dx=-1, maxiters=20, parallel= false, verb=0, atol=1e-10, n=6) =
                                    InterpolationMethod(dx, maxiters, parallel, verb, atol, n)

function findzero_interp(f::Function, x0::Float64, m::InterpolationMethod)
    dx = m.dx < 0 ? abs(x0)/20 : m.dx
    n = m.n
    s = x0
    ok = false
    iter = 1
    normf0 = Inf
    while !ok && iter <= m.maxiters
        m.verb > 3 && println("# TRIAL $iter for findzero_interp")
        xmax = s + 0.5*(n-1)*dx
        xmin  = s - 0.5*(n-1)*dx
        r = collect(linspace(xmin, xmax, n))
        if m.parallel
            refs = RemoteRef[]
            for i=1:length(r)
                push!(refs, @spawn f(r[i]))
            end

            f0 = [fetch(refs[i]) for i=1:length(r)]
        else
            f0 = [f(r[i]) for i=1:length(r)]
        end
        dummyfile, rf = mktemp()
        for i=1:length(r)
            println(rf, "$(r[i]) $(f0[i])")
        end
        close(rf)
        try
            s = float(readstring(`gnuplot -e "filename='$dummyfile'" -e "d=$s" $gplfile3`))
        catch
            error("ERROR GNUPLOT")
        end
        rm(dummyfile)

        normf0 = abs(f(s))
        if normf0 < m.atol
            m.verb > 3 && println("# SUCCESS x* = $(s), normf = $normf0")
            ok =true
        else
            m.verb > 1 && @warn("failed: iter=$iter x* = $(s), normf = $normf0")
            ok =false
        end

        if !ok && f0[1]*f0[end] > 0
            dx *= 2
            m.verb > 3 && println("# dx=$dx")
        elseif !ok && f0[1]*f0[end] < 0
            dx /= 2
            m.verb > 3 && println("# dx=$dx")
        end
        iter += 1
    end
    return ok, s, iter, normf0
end

max_stars_repo_name: JuliaTagBot/LittleScienceTools

max_issues_repo_name: JuliaTagBot/LittleScienceTools

max_forks_repo_name: JuliaTagBot/LittleScienceTools


text:
module StatisticalMachineLearning

using Statistics
using Random
import NearestNeighbors: knn
import Base: show

export Picture, TrainTestSplit, remove_constant, flatten, unflatten, visualize_picture, classify, knn, knn_acc, datamat

"""
    Picture(ID::Int, class::Int, data::Vector{T})

A type to hold the information about a picture.
"""
struct Picture{T}
    ID::Int
    class::Int
    data::AbstractVector{T}
    function Picture(x::AbstractVector{<:Real})
        try
            Int64(sqrt(length(x[begin+2:end])))
        catch e
            error("The input vector, with the first two values removed, is not square")
        end
        return new{eltype(x)}(x[1], x[2], x[3:end])
    end
    Picture(a::Int, b::Int, c::Vector{T}) where {T<:Real} = new{T}(a, b, c)
end

"""
    TrainTestSplit(pics::AbstractVector{Picture}, train_to_test_ratio::Rational, shuffle_pics::Bool=true)
    TrainTestSplit(ratio::Rational, n::Int, train::AbstractVector{Picture}, test::AbstractVector{Picture{T}})
"""
struct TrainTestSplit{T}
    ratio::Rational
    n::Int
    train::AbstractVector{Picture{T}}
    test::AbstractVector{Picture{T}}
    TrainTestSplit(ratio::Rational, n::Int, train::AbstractVector{Picture{T}}, test::AbstractVector{Picture{T}}) where {T<:Real} = new{T}(ratio, n, train, test)
    TrainTestSplit(train::AbstractVector{Picture{T}}, test::AbstractVector{Picture{T}}) where {T<:Real} = new{T}(length(train)//length(test), length(train)+length(test), train, test)
    function TrainTestSplit(pics::AbstractVector{Picture{T}}, train_to_test_ratio::Rational, shuffle_pics::Bool=true) where {T<:Real}
    if shuffle_pics
        shuffled_inds = shuffle(pics|>eachindex)
        pics = pics[shuffled_inds]
    end
    n_parts = train_to_test_ratio.num + train_to_test_ratio.den

    trainpics = pics[begin : train_to_test_ratio.num//n_parts * end ÷1]
    testpics = pics[(end - train_to_test_ratio.den//n_parts * end ÷1 +1) : end]
    return new{T}(train_to_test_ratio, length(pics), trainpics, testpics)
    end
end
testclasses(tts::TrainTestSplit) = getfield.(tts.test, :class)
trainclasses(tts::TrainTestSplit) = getfield.(tts.train, :class)

show(io::IO, p::Picture) = println(io, "A $(p.class) drawn by $(p.ID)")
show(io::IO, tts::TrainTestSplit) = println(io, "A TrainTestSplit object with $(tts.n) entries, of train to test ratio $(tts.ratio)")

"""
    flatten(m::AbstractMatrix) = vcat(eachcol(m)...)
"""
flatten(m::AbstractMatrix) = vcat(eachcol(m)...)
function unflatten(x::AbstractArray)
    try
        sidelength = sqrt(length(x)) |> Int64
    catch e
        "The sidelength (`x|>length|>sqrt`) could not be converted to an integer. "
    end
    batches = [x[sidelength*(i-1)+1:sidelength*i] for i in 1:sidelength]
    return hcat(batches...)
end

"""
    visualize_picture(data::AbstractMatrix, colormap::Symbol=:viridis)
    visualize_picture(data::AbstractVector, colormap::Symbol=:viridis)
    visualize_picture(p::Picture, colormap::Symbol=:viridis)
"""
function visualize_picture(data::AbstractMatrix, colormap::Symbol=:viridis)
    fig, ax = heatmap(data, figure=(resolution=(400, 400),); colormap)
    hidedecorations!(ax)
    fig
end
function visualize_picture(data::AbstractVector, colormap::Symbol=:viridis)
    fig, ax = heatmap(unflatten(data)|>reverse|>x->reverse(x, dims=1), figure=(resolution=(400, 400),); colormap)
    hidedecorations!(ax)
    fig
end
function visualize_picture(p::Picture, colormap::Symbol=:viridis)
    fig, ax = heatmap(unflatten(p.data)|>reverse|>x->reverse(x, dims=1), figure=(resolution=(400, 400),),
        axis=(title="Student ID: $(p.ID|>Int64)\nGround truth: $(p.class|>Int64)", ); colormap
    )
    hidedecorations!(ax)
    fig
end

"""
    classify(neighbor_inds::Vector{Int}, train_classes::Vector{Int}; tiebreaker=rand, possible_classes=unique(train_classes))

kwargs:
`tiebreaker` is
1) a function that takes a tuple of candidates and returns a value, or
2) a value that is returned upon a tie.
"""
function classify(neighbor_inds::Vector{Int}, train_classes::Vector{Int}, possible_classes=unique(train_classes); tiebreaker=rand)
    neighbor_classes = train_classes[neighbor_inds]
    my_counts = [count(==(possible_classes[i]), neighbor_classes) for i in eachindex(possible_classes)]
    A = [possible_classes my_counts]
    sorted_counts = sortslices(A, dims=1, by=x->x[2], rev=true)
    if sorted_counts[1, 2] == sorted_counts[2, 2]
        inds = [sorted_counts[i, 2] == sorted_counts[1, 2] for i in 1:size(sorted_counts, 1)]
        candidates_of_equal_count = sorted_counts[inds, :][:, 1]
        if tiebreaker isa Function
            return candidates_of_equal_count |> tiebreaker
        else
            f = (args...)->tiebreaker
            return candidates_of_equal_count |> f
        end
    else
        return sorted_counts[1, 1]
    end
end
function classify(neighbor_inds::Vector{Vector{Int}}, args...; kwargs...)
    [classify(neighbor_inds[i], args...; kwargs...) for i in eachindex(neighbor_inds)]
end

function knn(train_pics::Vector{Picture{T}}, test_pics::Vector{Picture{T}}; k::Int, tree=BruteTree, metric=Euclidean(), leafsize::Int=10) where {T<:Real}
    if tree == BruteTree
        mytree = tree(hcat(getfield.(train_pics, :data)...), metric)
    else
        mytree = tree(hcat(getfield.(train_pics, :data)...), metric; leafsize)
    end
    return knn(mytree, hcat(getfield.(test_pics, :data)...), k, true)
end
knn(tts::TrainTestSplit{<:Real}; kwargs...) = knn(tts.train, tts.test; kwargs...)

"""
    knn_acc(tts::TrainTestSplit{<:Real}; tiebreaker=rand, kwargs...)

`kwargs...` are passed on as `k::Int, tree=BruteTree, metric=Euclidean()`
"""
function knn_acc(tts::TrainTestSplit{<:Real}; tiebreaker=rand, kwargs...)
	inds, _ = knn(tts.train, tts.test; kwargs...)
	preds = classify(inds, trainclasses(tts); tiebreaker)
	return mean(preds .== testclasses(tts))
end

"""
	knn_acc_crossvalidate(pics::Vector{Picture}, ratio::Rational = 9//1)

Shuffle `pics` ratio.num + ratio.den times, make a TrainTestSplit for each, and calculate
the mean and standard deviation of the classification accuracy of each.
"""
function knn_acc_crossvalidate(pics::Vector{Picture{T}}, ratio::Rational = 9//1) where {T<:Real}
	n_runs = ratio.num + ratio.den
	temp_TTSs = [TrainTestSplit(pics, ratio) for _ in 1:n_runs]
    x = [knn_acc(tts, k=3) for tts in temp_TTSs]
    (mean=mean(x), std=std(x), n_runs = n_runs)
end

"""
    datamat(pics::Vector{Picture}) = hcat(getfield.(pics, :data)...)

Return the data of all pictures in `pics`.
A column in the returned Matrix represents a single picture.
"""
datamat(pics::Vector{Picture{T}}) where {T<:Real} = hcat(getfield.(pics, :data)...)
datamat(tts::TrainTestSplit{T}) where {T<:Real} = (train=datamat(tts.train), test=datamat(tts.test))

function remove_constant(m::AbstractMatrix)
    bad_row_inds = Int64[]
    all_row_inds = 1:size(m, 1) |> Vector
    for i in all_row_inds
        if maximum(m[i, :]) == minimum(m[i, :])
            push!(bad_row_inds, i)
        end
    end
    reduced_data = m[deleteat!(all_row_inds, bad_row_inds), :]
    return reduced_data
end

function remove_constant(pics::Vector{Picture{T}}) where {T<:Real}
    datatogether = datamat(pics)
    reduced_datamat = remove_constant(datatogether)
    return [Picture(pics[i].ID, pics[i].class, reduced_datamat[:, i]) for i in eachindex(pics)]
end

end


text:
"""
    filter_non_term_leaves!(g, term_vert)

Remove edges of `g` so that all non-isolated leaves of `g` are in the set `term_vert`
"""
function filter_non_term_leaves!(
    g::AbstractGraph{T},
    term_vert::Vector{<:Integer}
    ) where T<:Integer

    is_term = falses(nv(g))
    is_term[term_vert] .= true
    leaves = [v for v in vertices(g) if degree(g, v) == 1 && !is_term[v]]

    while !isempty(leaves)
        leaf = pop!(leaves)
        for v in neighbors(g, leaf)
            rem_edge!(g, v, leaf)
            if degree(g, v) == 1 && !is_term[v]
                push!(leaves, v)
            end
        end
    end

    return g
end

"""
	steiner_tree(g, term_vert, distmx=weights(g))

Return an approximately minimum steiner tree of connected, undirected graph `g` with positive edge
weights represented by `distmx` using [Approximate Steiner Tree](https://en.wikipedia.org/wiki/Steiner_tree_problem#Approximating_the_Steiner_tree).
The minimum steiner tree problem involves finding a subset of edges in `g` of minimum weight such
that all the vertices in `term_vert` are connected.

`t = length(term_vert)`.

### Performance
Runtime: O(t*(t*log(t)+|E|*log(|V| ))
Memory: O(t*|V|)
Approximation Factor: 2-2/t
"""
function steiner_tree end

@traitfn function steiner_tree(
    g::AG::(!IsDirected),
    term_vert::Vector{<:Integer},
    distmx::AbstractMatrix{U} = weights(g)
    ) where {U<:Real, T, AG<:AbstractGraph{T}}

    nvg = nv(g)
    term_to_actual = T.(term_vert)
    unique!(term_to_actual)

    # Compute the graph formed by inducing the metric closure graph of g on term_vert
    nvg_mc = length(term_to_actual)
    distmx_mc = Matrix{U}(undef, nvg_mc, nvg_mc)
    parents = Matrix{T}(undef, nvg, nvg_mc)

    for (i, v) in enumerate(term_to_actual)
        d_s = dijkstra_shortest_paths(g, v, distmx)
        @inbounds parents[:, i] = d_s.parents
        @inbounds distmx_mc[:, i] = @view d_s.dists[term_to_actual]
    end

    # Compute MST of Metric Closure graph
    mst_mc = kruskal_mst(CompleteGraph(nvg_mc), distmx_mc)
    expanded_mst = Vector{Edge{T}}()
    sizehint!(expanded_mst, nvg)

    # Expand each edge in mst_mc into the path it represents
    for e in mst_mc
        i = src(e)
        s = term_to_actual[i]
        t = term_to_actual[dst(e)]
        while s != t
            t_next = parents[t, i]
            push!(expanded_mst, Edge(min(t_next, t), max(t_next, t)))
            t = t_next
        end
    end

    # Compute the MST of the expanded graph
    mst_mst_mc = kruskal_mst(SimpleGraph(expanded_mst), distmx)

    # Remove non-terminal leaves
    return filter_non_term_leaves!(SimpleGraph(mst_mst_mc), term_to_actual)
end


text:
include("../../../experiments/OceanSplitExplicit/simple_box.jl")
ClimateMachine.init()

# Float type
const FT = Float64

#################
# RUN THE TESTS #
#################
@testset "$(@__FILE__)" begin

    include("../refvals/simple_box_ivd_refvals.jl")
    refDat = (refVals[1], refPrecs[1])

    # simulation time
    timestart = FT(0)      # s
    timeend = FT(5 * 86400) # s
    timespan = (timestart, timeend)

    # DG polynomial order
    N = Int(4)

    # Domain resolution
    Nˣ = Int(20)
    Nʸ = Int(20)
    Nᶻ = Int(20)
    resolution = (N, Nˣ, Nʸ, Nᶻ)

    # Domain size
    Lˣ = 4e6    # m
    Lʸ = 4e6    # m
    H = 1000   # m
    dimensions = (Lˣ, Lʸ, H)

    BC = (
        ClimateMachine.Ocean.SplitExplicit01.CoastlineNoSlip(),
        ClimateMachine.Ocean.SplitExplicit01.OceanFloorNoSlip(),
        ClimateMachine.Ocean.SplitExplicit01.OceanSurfaceStressForcing(),
    )

    config = config_simple_box(
        "test_simple_box",
        resolution,
        dimensions,
        BC;
        dt_slow = FT(90 * 60),
        dt_fast = FT(240),
    )

    run_simple_box(config, timespan; refDat = refDat)

end

max_stars_repo_name: mwarusz/CLIMA

max_issues_repo_name: mwarusz/CLIMA

max_forks_repo_name: mwarusz/CLIMA


text:
function TestJacobi()
dlambda=pi/8;
dtheta=pi/8;
s1=[0 0];
s2=[dlambda 0];
s3=[dlambda dtheta];
s4=[0 dtheta];

c1=SphereToCart(s1(1),s1(2));
c2=SphereToCart(s2(1),s2(2));
c3=SphereToCart(s3(1),s3(2));
c4=SphereToCart(s4(1),s4(2));

xi1=.3;
xi2=.4;
cTilde=rTilde(xi1,xi2,c1,c2,c3,c4);
c=cTilde/norm(cTilde);
[lambda,theta]=CartToSphereToCart(c);

cNeu=SphereToCart(lambda,theta);

DV1=0.25/norm(cTilde)*JacS1(lambda,theta)*JacS2(lambda,theta);

D0=[1/norm(cTilde) 0
   0 norm(cTilde)^2];

 DV2=0.25*D0*JacC1(cTilde)*JacC2(cTilde);

 DV3=0.25/norm(cTilde)^(3/2)*JacCC1(c)*JacCC2(cTilde);
 ee=dd();
end

function x=SphereToCart(lambda,theta)
x=zeros(3,1);
x(1)=cos(lambda)*cos(theta);
x(2)=sin(lambda)*cos(theta);
x(3)=sin(theta);
end

function [lambda,theta]=CartToSphereToCart(x)
r=norm(x);
lambda=atan2(x(2),x(1));
if x(2)<0
  lambda=lambda+2*pi;
end
theta=asin(x(3)/r);
end

function x=rTilde(xi1,xi2,c1,c2,c3,c4)
x=0.25*((1-xi1)*(1-xi2)*c1...
  +(1+xi1)*(1-xi2)*c2...
  +(1+xi1)*(1+xi2)*c3...
  +(1-xi1)*(1+xi2)*c4);
end

function D=JacCC1(x)
D=zeros(2,2);
D(1,1)=-x(2)/(x(1)*x(1)+x(2)*x(2));
D(1,2)=x(1)/(x(1)*x(1)+x(2)*x(2));

D(2,1)=x(1)/((x(1)^2 + x(2)^2)^(1/2)*(- x(1)^2 - x(2)^2 + 1)^(1/2));
D(2,2)=x(2)/((x(1)^2 + x(2)^2)^(1/2)*(- x(1)^2 - x(2)^2 + 1)^(1/2));
end

function D=JacCC2(x)
D=zeros(2,3);
D(1,1)=x(2)*x(2)+x(3)*x(3);
D(1,2)=-x(1)*x(2);
D(1,3)=-x(1)*x(3);

D(2,2)=-x(1)*x(2);
D(2,2)=x(1)*x(1)+x(3)*x(3);
D(2,3)=-x(2)*x(3);

end

function D=JacC1(x)
D=zeros(2,3);
D(1,1)=-x(2)/(x(1)*x(1)+x(2)*x(2));
D(1,2)=x(1)/(x(1)*x(1)+x(2)*x(2));

D(2,1)=x(1)*x(3)/sqrt(x(1)*x(1)+x(2)*x(2));
D(2,2)=x(2)*x(3)/sqrt(x(1)*x(1)+x(2)*x(2));
D(2,3)=-sqrt(x(1)*x(1)+x(2)*x(2));
end

function D=JacC2(x)
D=zeros(3,3);
D(1,1)=x(2)*x(2)+x(3)*x(3);
D(1,2)=-x(1)*x(2);
D(1,3)=-x(1)*x(3);

D(2,2)=-x(1)*x(2);
D(2,2)=x(1)*x(1)+x(3)*x(3);
D(2,3)=-x(2)*x(3);

D(3,1)=-x(1)*x(3);
D(3,2)=-x(2)*x(3);
D(3,3)=x(1)*x(1)+x(2)*x(2);
end

function D=JacS1(lambda,theta)
D=zeros(2,3);
D(1,1)=-sin(lambda);
D(1,2)=cos(lambda);

D(2,3)=1;
end

function D=JacS2(lambda,theta)
D=zeros(3,3);
D(1,1)=sin(lambda)^2*cos(theta)^2+sin(theta)^2;
D(1,2)=-0.5*sin(2*lambda)*cos(theta)^2;
D(1,3)=0.5*cos(lambda)*sin(2*theta);

D(2,1)=-0.5*sin(2*lambda)*cos(theta)^2;
D(2,2)=cos(lambda)^2*cos(theta)^2+sin(theta)^2;
D(2,3)=-0.5*sin(lambda)*sin(2*theta);

D(3,1)=-0.5*cos(lambda)*sin(theta);
D(3,2)=-0.5*sin(lambda)*sin(theta);
D(3,3)=cos(theta);
end

function erg=dd()
syms x
syms y
f=acos(sqrt(1-x*x-y*y));
erg=diff(f,x);
end


text:
using DiffEqOperators, Setfield, Parameters
using LinearAlgebra, Plots, SparseArrays
using BifurcationKit

plotsol(x, Nx=Nx, Ny=Ny) = heatmap(reshape(Array(x), Nx, Ny)', color=:viridis)
plotsol!(x, Nx=Nx, Ny=Ny; kwargs...) = heatmap!(reshape(Array(x), Nx, Ny)'; color=:viridis, kwargs...)

Nx = 151
	Ny = 100
	lx = 8pi
	ly = 2*2pi/sqrt(3)

function Laplacian2D(Nx, Ny, lx, ly, bc = :Neumann)
	hx = 2lx/Nx
	hy = 2ly/Ny
	D2x = CenteredDifference(2, 2, hx, Nx)
	D2y = CenteredDifference(2, 2, hy, Ny)
	if bc == :Dirichlet
		Qx = Dirichlet0BC(typeof(hx))
		Qy = Dirichlet0BC(typeof(hy))
	elseif bc == :Neumann
		Qx = Neumann0BC(hx)
		Qy = Neumann0BC(hy)
	elseif bc == :Periodic
		Qx = PeriodicBC(hx)
		Qy = PeriodicBC(hy)
	end
	# @show norm(D2x - Circulant(D2x[1,:]))
	A = kron(sparse(I, Ny, Ny), sparse(D2x * Qx)[1]) + kron(sparse(D2y * Qy)[1], sparse(I, Nx, Nx))
	return A, D2x
end

function F_sh(u, p)
	@unpack l, ν, L1 = p
	return -L1 * u .+ (l .* u .+ ν .* u.^2 .- u.^3)
end

function dF_sh(u, p)
	@unpack l, ν, L1 = p
	return -L1 .+ spdiagm(0 => l .+ 2 .* ν .* u .- 3 .* u.^2)
end

d2F_sh(u, p, dx1, dx2) = (2 .* p.ν .* dx2 .- 6 .* dx2 .* u) .* dx1
d3F_sh(u, p, dx1, dx2, dx3) = (-6 .* dx2 .* dx3) .* dx1
jet = (F_sh, dF_sh, d2F_sh, d3F_sh)

X = -lx .+ 2lx/(Nx) * collect(0:Nx-1)
Y = -ly .+ 2ly/(Ny) * collect(0:Ny-1)

sol0 = [(cos(x) .+ cos(x/2) * cos(sqrt(3) * y/2) ) for x in X, y in Y]
	sol0 .= sol0 .- minimum(vec(sol0))
	sol0 ./= maximum(vec(sol0))
	sol0 = sol0 .- 0.25
	sol0 .*= 1.7
	heatmap(sol0', color=:viridis)

Δ, D2x = Laplacian2D(Nx, Ny, lx, ly, :Neumann)
L1 = (I + Δ)^2
par = (l = -0.1, ν = 1.3, L1 = L1)

optnew = NewtonPar(verbose = true, tol = 1e-8, maxIter = 20)
# optnew = NewtonPar(verbose = true, tol = 1e-8, maxIter = 20, eigsolver = EigArpack(0.5, :LM))
	sol_hexa, hist, flag = @time newton(F_sh, dF_sh, vec(sol0), par, optnew)
	println("--> norm(sol) = ", norm(sol_hexa, Inf64))
	plotsol(sol_hexa)


text:
using HTTPClient.HTTPC
using JSON

api = "https://www.googleapis.com/books/v1/volumes"
url = api * "?q=isbn:1408855895"
uu = HTTPC.get(url)
json = JSON.parse(bytestring(uu.body.data));

volumeInfo = json["items"][1]["volumeInfo"]
title   = volumeInfo["title"]
author  = volumeInfo["authors"][1]
publisher = volumeInfo["publishedBy"]
pubdate = volumeInfo["publishedDate"]
ppcoount = volumeInfo["pageCount"]

@printf "%s by %s\nPublished by %s (%s), Pages: $d" title, author, publisher, pubdate, ppcount

max_stars_repo_path: Module 3/Chapter08/rest-books.jl

max_issues_repo_path: Module 3/Chapter08/rest-books.jl

max_forks_repo_path: Module 3/Chapter08/rest-books.jl


text:
"""
Parse Elements of Map
"""
function parse_element(handler::LibExpat.XPStreamHandler,
                      name::AbstractString,
                      attr::Dict{AbstractString,AbstractString})
    data = handler.data::OpenStreetMapX.DataHandle
    if name == "node"
        data.element = :Tuple
        data.node = (parse(Int, attr["id"]),
                         OpenStreetMapX.LLA(parse(Float64,attr["lat"]), parse(Float64,attr["lon"])))
    elseif name == "way"
        data.element = :Way
        data.way = OpenStreetMapX.Way(parse(Int, attr["id"]))
    elseif name == "relation"
        data.element = :Relation
        data.relation = OpenStreetMapX.Relation(parse(Int, attr["id"]))
    elseif name == "bounds"
        data.element =:Bounds
        data.bounds = OpenStreetMapX.Bounds(parse(Float64,attr["minlat"]), parse(Float64,attr["maxlat"]), parse(Float64,attr["minlon"]), parse(Float64,attr["maxlon"]))
    elseif name == "tag"
        k = attr["k"]; v = attr["v"]
        if  data.element == :Tuple
			if haskey(FEATURE_CLASSES, k)
				data.osm.features[handler.data.node[1]] = k,v
			end
		elseif data.element == :Way
            data_tags = tags(data.way)
            push!(data.osm.way_tags, k)
			data_tags[k] = v
        elseif data.element == :Relation
            data_tags = tags(data.relation)
            push!(data.osm.relation_tags, k)
			data_tags[k] = v
        end
    elseif name == "nd"
        push!(data.way.nodes, parse(Int, attr["ref"]))
    elseif name == "member"
        push!(data.relation.members, attr)
    end
end

function collect_element(handler::LibExpat.XPStreamHandler, name::AbstractString)
    if name == "node"
        handler.data.osm.nodes[handler.data.node[1]] = handler.data.node[2]
        handler.data.element = :None
    elseif name == "way"
        push!(handler.data.osm.ways, handler.data.way)
        handler.data.element = :None
    elseif name == "relation"
        push!(handler.data.osm.relations, handler.data.relation)
        handler.data.element = :None
    elseif name == "bounds"
        handler.data.osm.bounds = handler.data.bounds
		handler.data.element = :None
    end
end

function parseOSM(filename::AbstractString; args...)::OpenStreetMapX.OSMData
    callbacks = LibExpat.XPCallbacks()
    callbacks.start_element = parse_element
    callbacks.end_element = collect_element
    data = OpenStreetMapX.DataHandle()
    LibExpat.parsefile(filename, callbacks, data=data; args...)
    data.osm
end

"""
    get_map_data(filepath::String,filename::Union{String,Nothing}=nothing;
	             road_levels::Set{Int} = Set(1:length(OpenStreetMapX.ROAD_CLASSES)),
				 use_cache::Bool = true, only_intersections::Bool=true)::MapData

High level function - parses .osm file and create the road network based on the map data.

**Arguments**

* `filepath` : path with an .osm file (directory or path to a file)
* `filename` : name of the file (when the first argument is a directory)
* `road_levels` : a set with the road categories (see: OpenStreetMapX.ROAD_CLASSES for more informations)
* `use_cache` : a *.cache file will be crated with a serialized map image in the `datapath` folder
* `only_intersections` : include only road system data
* `trim_to_connected_graph`: trim orphan nodes in such way that the map is a strongly connected graph
"""
function get_map_data(filepath::String,filename::Union{String,Nothing}=nothing; road_levels::Set{Int} = Set(1:length(OpenStreetMapX.ROAD_CLASSES)),
		use_cache::Bool = true,only_intersections::Bool = true, trim_to_connected_graph::Bool=false)::MapData
    #preprocessing map file
    datapath = (filename==nothing) ? dirname(filepath) : filepath;
	if filename == nothing
		filename = basename(filepath)
	end
	cachefile = joinpath(datapath,filename*".cache")
	if use_cache && isfile(cachefile)
		f=open(cachefile,"r");
		res=Serialization.deserialize(f);
		close(f);
		@info "Read map data from cache $cachefile"
	else
		mapdata = OpenStreetMapX.parseOSM(joinpath(datapath,filename))
		OpenStreetMapX.crop!(mapdata,crop_relations = false)
		res = MapData(mapdata, road_levels, only_intersections; trim_to_connected_graph=trim_to_connected_graph)
		if use_cache
			f=open(cachefile,"w");
			Serialization.serialize(f,res);
			@info "Saved map data to cache $cachefile"
			close(f);
		end
	end
    return res
end

"""
Get Vertices and nodes for a set of `edges`
"""
function get_vertices_and_graph_nodes(edges::Vector{Tuple{Int,Int}})
    graph_nodes = unique(reinterpret(Int, edges))
	vertices = Dict{Int,Int}(zip(graph_nodes, 1:length(graph_nodes)))
	return vertices, graph_nodes
end

"""
Internal constructor of `MapData` object
"""
function MapData(mapdata::OSMData, road_levels::Set{Int}, only_intersections::Bool=true;
	 trim_to_connected_graph::Bool=false, remove_nodes::AbstractSet{Int}=Set{Int}())
	#preparing data
	bounds = mapdata.bounds
	nodes = OpenStreetMapX.ENU(mapdata.nodes,OpenStreetMapX.center(bounds))
	highways = OpenStreetMapX.filter_highways(OpenStreetMapX.extract_highways(mapdata.ways))
	roadways = OpenStreetMapX.filter_roadways(highways, levels= road_levels)
	if length(remove_nodes) > 0
		delete!.(Ref(nodes), remove_nodes);
		delcount = 0
		for rno in length(roadways):-1:1
			rr = roadways[rno]
			for i in length(rr.nodes):-1:1
				if rr.nodes[i] in remove_nodes
					deleteat!(rr.nodes,i)
					delcount += 1
				end
			end
			length(rr.nodes) == 0 && deleteat!(roadways, rno)
		end
	end
	intersections = OpenStreetMapX.find_intersections(roadways)
	segments = OpenStreetMapX.find_segments(nodes,roadways,intersections)
	#remove unuseful nodes
	roadways_nodes = unique(vcat(collect(way.nodes for way in roadways)...))
	nodes = Dict(key => nodes[key] for key in roadways_nodes)

	# e - Edges in graph, stored as a tuple (source,destination)
	# class - Road class of each edgey
	if only_intersections && !trim_to_connected_graph
		vals = Dict((segment.node0,segment.node1) => (segment.distance,segment.parent) for segment in segments)
		e = collect(keys(vals))
		vals = collect(values(vals))
		weight_vals = map(val -> val[1],vals)
		classified_roadways = OpenStreetMapX.classify_roadways(roadways)
		class =  [classified_roadways[id] for id in map(val -> val[2],vals)]
	else
		e,class = OpenStreetMapX.get_edges(nodes,roadways)
		weight_vals = OpenStreetMapX.distance(nodes,e)
	end
	# (node id) => (graph vertex)
	v, n = OpenStreetMapX.get_vertices_and_graph_nodes(e)
	edges = [v[id] for id in reinterpret(Int, e)]
	I = edges[1:2:end]
	J = edges[2:2:end]
	# w - Edge weights, indexed by graph id
	w = SparseArrays.sparse(I, J, weight_vals, length(v), length(v))
	g = LightGraphs.DiGraph(length(v))
	for edge in e
		add_edge!(g,v[edge[1]], v[edge[2]])
	end

	if trim_to_connected_graph
		rm_list = Set{Int}()
		conn_components = sort!(LightGraphs.strongly_connected_components(g),
        	lt=(x,y)->length(x)<length(y), rev=true)
		remove_vs = vcat(conn_components[2:end]...)
		rm_list = getindex.(Ref(n), remove_vs)
		return  MapData(mapdata, road_levels, only_intersections, remove_nodes=Set{Int}(rm_list))
	else
		return MapData(bounds,nodes,roadways,intersections,g,v,n,e,w,class)
	end
end


text:
# sort by column
function sortbyc(x,col)
    x = sortslices(x,dims=1,by=x->x[col])
end

max_stars_repo_name: nilshg/Econometrics

max_issues_repo_name: nilshg/Econometrics

max_forks_repo_name: nilshg/Econometrics


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

struct ApisApi <: SwaggerApi
    client::Swagger.Client
end

"""

get available API versions
Return: IoK8sApimachineryPkgApisMetaV1APIGroupList
"""
function getAPIVersions(_api::ApisApi; _mediaType=nothing)
    _ctx = Swagger.Ctx(_api.client, "GET", IoK8sApimachineryPkgApisMetaV1APIGroupList, "/apis/", ["BearerToken"])
    Swagger.set_header_accept(_ctx, ["application/json", "application/yaml", "application/vnd.kubernetes.protobuf"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? ["application/json", "application/yaml", "application/vnd.kubernetes.protobuf"] : [_mediaType])
    Swagger.exec(_ctx)
end

export getAPIVersions


text:
"""
    ProviderEHRLaunchConfig(; kwargs...)

## Required Keyword Arguments:
- `client_id::String`
- `redirect_uri::String`
"""
Base.@kwdef struct ProviderEHRLaunchConfig
    client_id::String
    redirect_uri::String
end

"""
    ProviderStandaloneConfig(; kwargs...)

## Required Keyword Arguments:
- `authorize_endpoint::String`
- `client_id::String`
- `redirect_uri::String`
- `token_endpoint::String`
"""
Base.@kwdef struct ProviderStandaloneConfig
    authorize_endpoint::String
    client_id::String
    redirect_uri::String
    token_endpoint::String
end

Base.@kwdef struct _AuthorizationCodeInformation
    authorization_code::String
    client_id::String
    redirect_uri::String
    token_endpoint::String
end

Base.@kwdef struct _AccessTokenInformation
    access_token::String
    access_token_is_jwt::Bool = false
    access_token_jwt_decoded::Dict{String, Any} = Dict{String, Any}()
    access_token_response::Dict{Symbol, Any}
    authorization_code::String
    authorization_code_is_jwt::Bool = false
    authorization_code_jwt_decoded::Dict{String, Any} = Dict{String, Any}()
end

Base.@kwdef struct ProviderEHRLaunchResult
    access_token::String
    access_token_is_jwt::Bool = false
    access_token_jwt_decoded::Dict{String, Any} = Dict{String, Any}()
    access_token_response::Dict{Symbol, Any}
    authorization_code::String
    authorization_code_is_jwt::Bool = false
    authorization_code_jwt_decoded::Dict{String, Any} = Dict{String, Any}()
    launch_token::String
    launch_token_is_jwt::Bool = false
    launch_token_jwt_decoded::Dict{String, Any} = Dict{String, Any}()
end

Base.@kwdef struct ProviderStandaloneResult
    access_token::String
    access_token_is_jwt::Bool = false
    access_token_jwt_decoded::Dict{String, Any} = Dict{String, Any}()
    access_token_response::Dict{Symbol, Any}
    authorization_code::String
    authorization_code_is_jwt::Bool = false
    authorization_code_jwt_decoded::Dict{String, Any} = Dict{String, Any}()
end

"""
    HealthBase.get_fhir_access_token(result::ProviderEHRLaunchResult)
"""
function HealthBase.get_fhir_access_token(result::ProviderEHRLaunchResult)
    return result.access_token
end

"""
    HealthBase.has_fhir_patient_id(result::ProviderEHRLaunchResult)
"""
function HealthBase.has_fhir_patient_id(result::ProviderEHRLaunchResult)
    return haskey(result.access_token_response, :patient)
end

"""
    HealthBase.get_fhir_patient_id(result::ProviderEHRLaunchResult)
"""
function HealthBase.get_fhir_patient_id(result::ProviderEHRLaunchResult)
    return result.access_token_response[:patient]
end

"""
    HealthBase.get_fhir_access_token(result::ProviderStandaloneResult)
"""
function HealthBase.get_fhir_access_token(result::ProviderStandaloneResult)
    return result.access_token
end


text:
module day13

import ..aoc19: run_program!
using StaticArrays, Combinatorics, ProgressMeter

readInput() = parse.(Int, split(readline(joinpath(@__DIR__,"input.txt")), ','))

struct Out
    x::Int
    y::Int
    id::Int
end

function part1(data)
    input = Channel{Int}(Inf)
    output = Channel{Int}(Inf)
    run_program!(data, input, output)
    close(output)
    result = [o for o in output]
    count(==(2), (result[i] for i in 3:3:length(result)))
end

function part2(data)
    data[1] = 2 # insert 2 quarters

    input = Channel{Int}(Inf)
    output = Channel{Int}(Inf)
    game = @async run_program!(data, input, output)

    ball = Out(0, 0, 0)
    paddle = Out(0, 0, 0)
    score = Out(0, 0, 0)

    while !istaskdone(game)
        yield()
        result = Out[]
        while isready(output)
            push!(result, Out(take!(output), take!(output), take!(output)))
        end

        for out in result
            if out.id == 4
                ball = out
            elseif out.id == 3
                paddle = out
            elseif out.x == -1 && out.y == 0
                score = out
            end
        end

        if ball.x < paddle.x
            put!(input, -1)
        elseif ball.x > paddle.x
            put!(input, 1)
        else
            put!(input, 0)
        end
    end

    score.id
end

end


text:
"""
    AbstractEnclosureAlgorithm

Abstract type for range enclosure algorithms.
"""
abstract type AbstractEnclosureAlgorithm end

"""
   NaturalEnclosure <: AbstractEnclosureAlgorithm

Data type to bound the range of `f` over `X` using natural enclosure, i.e., to
evaluate `f(X)` with interval arithmetic.

### Examples

```jldoctest
julia> enclose(x -> 1 - x^4 + x^5, 0..1, NaturalEnclosure())
[0, 2]
```
"""
struct NaturalEnclosure <: AbstractEnclosureAlgorithm end

"""
    MooreSkelboeEnclosure{T} <: AbstractEnclosureAlgorithm

Data type to bound the range of `f` over `X` using the Moore-Skelboe algorithm, which
rigorously computes the global minimum and maximum of the function.
See [`IntervalOptimisation.jl`](https://github.com/JuliaIntervals/IntervalOptimisation.jl) for more details.

### Fields

- `structure` -- (default: `HeapedVector`) the way in which vector elements
                 are kept arranged; possible options are `HeapedVector` and `SortedVector`
- `tol`       -- (default `1e-3`) tolerance to which the optima are computed

### Notes

To use this algorithm, you need to import `IntervalOptimisation.jl`

### Examples

```jldoctest
julia> using IntervalOptimisation

julia> enclose(x -> 1 - x^4 + x^5, 0..1, MooreSkelboeEnclosure()) # default parameters
[0.916034, 1.00213]

julia> enclose(x -> 1 - x^4 + x^5, 0..1, MooreSkelboeEnclosure(; tol=1e-2))
[0.900812, 1.0326]
"""
Base.@kwdef struct MooreSkelboeEnclosure{T} <: AbstractEnclosureAlgorithm
    structure::T = HeapedVector
    tol::Float64 = 1e-3
end

"""
    TaylorModelsEnclosure <: AbstractEnclosureAlgorithm

Data type to bound the range of `f` over `X` using Taylor models.
See [`TaylorModels.jl`](https://github.com/JuliaIntervals/TaylorModels.jl) for more details.

### Fields

- `order`     -- (default: `10`) order of the Taylor model used to compute
                 an enclosure of `f` over `dom`
- `normalize` -- (default: `true`) if `true`, normalize the Taylor model
                 on the unit symmetric box around the origin

### Examples

```jldoctest
julia> enclose(x -> 1 - x^4 + x^5, 0..1, TaylorModelsEnclosure()) # default parameters
[0.8125, 1.09375]

julia> enclose(x -> 1 - x^4 + x^5, 0..1, TaylorModelsEnclosure(; order=4))
[0.78125, 1.125]

```
"""
Base.@kwdef struct TaylorModelsEnclosure <: AbstractEnclosureAlgorithm
    order::Int = 10
    normalize::Bool = true
end

"""
    SumOfSquaresEnclosure{T} <: AbstractEnclosureAlgorithm

Data type to bound the range of `f` over `X` using sum-of-squares optimization.
See [`SumOfSquares.jl`](https://github.com/jump-dev/SumOfSquares.jl) for more details

### Fields

- `backend` -- backend used to solve the optimization problem; a list of available backends
             can be found [here](https://jump.dev/JuMP.jl/stable/installation/#Supported-solvers)
- `order` -- (default `5`), maximum degree of the SDP relaxation

### Notes

To use this solver, you need to import `SumOfSquares.jl` and a backend.
Since the optimization problem is solved numerically and not with interval arithmetic, the
result of this algorithm is not rigorous.

```julia
julia> using SumOfSquares, SDPA

julia> backend = SDPA.Optimizer

julia> enclose(x -> -x^3/6 + 5x, 1..4, SumOfSquaresEnclosure(; backend=backend))
[4.8333, 10.541]
```
"""
Base.@kwdef struct SumOfSquaresEnclosure{T} <: AbstractEnclosureAlgorithm
    backend::T
    order::Int = 5
end

"""
    BranchAndBoundEnclosure

Data type to bound the range of `f` over `X` using the branch and bound algorithm.

### Fields

- `maxdepth` (default `10`): maximum depth of the search tree
- `tol` (default `1e-3`): tolerance to compute the range of the function

### Algorithm

The algorithm evaluates a function `f` over an interval `X`. If the maximum depth is reached or the
width of `f(X)` is below the tolerance, the algorithm returns the computed range; otherwise it bisects
the interval/interval box.

The algorithm also looks at the sign of the derivative / gradient to see if the range can be
computed directly. By default, the derivative / gradient is computed using `ForwardDiff.jl`,
but a custom value can be passed via the `df` keyword argument to [`enclose`](@ref).

### Examples

```jldoctest
julia> enclose(x -> -x^3/6 + 5x, 1..4, BranchAndBoundEnclosure())
[4.83333, 10.5709]

julia> enclose(x -> -x^3/6 + 5x, 1..4, BranchAndBoundEnclosure(tol=1e-2); df=x->-x^2/2+5)
[4.83333, 10.5709]
"""
Base.@kwdef struct BranchAndBoundEnclosure <: AbstractEnclosureAlgorithm
    maxdepth = 10
    tol = 1e-3
end