text:
using Symbolics
using RoboDojo
using LinearAlgebra
using SparseArrays

using MeshCat
using GeometryBasics
using Colors
using CoordinateTransformations
using Rotations

# visualize
function set_background!(vis::Visualizer; top_color=RGBA(1,1,1.0), bottom_color=RGBA(1,1,1.0))
    setprop!(vis["/Background"], "top_color", top_color)
    setprop!(vis["/Background"], "bottom_color", bottom_color)
end
vis = Visualizer()
open(vis)
set_background!(vis)

function hat(x)
    return [0 -x[3] x[2];
            x[3] 0 -x[1];
           -x[2] x[1] 0]
end

function L_mult(x)
    [x[1] -transpose(x[2:4]);
     x[2:4] x[1] * I(3) + hat(x[2:4])]
end

# right quaternion multiply as matrix
function R_mult(x)
    [x[1] -transpose(x[2:4]); x[2:4] x[1] * I(3) - hat(x[2:4])]
end

# rotation matrix
function rotation_matrix(q)
    H = [zeros(1, 3); I(3)]
    transpose(H) * L_mult(q) * transpose(R_mult(q)) * H
end
"""
    variables:
    pa ∈ R^3 # parent body contact point
    pb ∈ R^3 # child  body  contact point

    sa ∈ R+  # slack for parent body collision set
    sb ∈ R+  # slack for child  body collision set

    λa ∈ R   # dual  for parent equality constraint
    λb ∈ R   # dual  for child  equality constraint

    νa ∈ R+  # dual  for parent slack
    νb ∈ R+  # dual  for child  slack

    data:
    xa ∈ R^3 # parent body position
    qa ∈ H   # parent body orientation
    xb ∈ R^3 # child  body position
    qb ∈ H   # child  body orientation
    ra ∈ R+  # parent body collision radius
    rb ∈ R+  # child  body collision radius

    objective:
    ||pa - pb||

    -> (pa - pb)' (pa - pb) # squared distance

    constraints:
    ||R(qa)' * (pa - xa)||_Σa^-1 <= 1.0 # ellipsoid
    ||R(qb)' * (pb - xb)||_Σb^-1 <= 1.0 # ellipsoid
        where Σ^-1 = Diagonal([1.0 / x^2, 1.0 / y^2, 1.0 / z^2])
    ->
        sa - (1.0 - (pa - xa)' R(qa) Σa^-1 R(qa)' (pa - xa)) = 0 # squared + slack
        sb - (1.0 - (pb - xb)' R(qb) Σb^-1 R(qb)' (pb - xb)) = 0 # squared + slack
        sa, sb >= 0

    primal:
    z = (pa, pb, sa, sb) ∈ R^8
    dual:
    y = (λa, λb, νa, νb) ∈ R^4
    data:
    θ = (xa, qa, xb, qb, ax, ay, az, bx, by, bz) ∈ R^16
"""

nz = 8
ny = 4
nw = nz + ny
nθ = 20
@variables z[1:nz] y[1:ny] θ[1:nθ] κ[1:1]

function unpack_primals(z)
    pa = z[0 .+ (1:3)]
    pb = z[3 .+ (1:3)]
    sa = z[6 .+ (1:1)]
    sb = z[7 .+ (1:1)]

    return pa, pb, sa, sb
end

function unpack_duals(y)
   λa = y[0 .+ (1:1)]
   λb = y[1 .+ (1:1)]
   νa = y[2 .+ (1:1)]
   νb = y[3 .+ (1:1)]

   return λa, λb, νa, νb
end

function unpack_data(θ)
    xa = θ[1:3]
    qa = θ[3 .+ (1:4)]
    xb = θ[7 .+ (1:3)]
    qb = θ[10 .+ (1:4)]
    ax = θ[14 .+ (1:1)]
    ay = θ[15 .+ (1:1)]
    az = θ[16 .+ (1:1)]
    bx = θ[17 .+ (1:1)]
    by = θ[18 .+ (1:1)]
    bz = θ[19 .+ (1:1)]

    return xa, qa, xb, qb, ax, ay, az, bx, by, bz
end

function objective(z)
    pa, pb, sa, sb = unpack_primals(z)
    d = pa - pb
    return dot(d, d)
end

function constraints(z, θ)
    pa, pb, sa, sb = unpack_primals(z)
    xa, qa, xb, qb, ax, ay, az, bx, by, bz = unpack_data(θ)
    [
        sa[1] - (1.0^2 - (pa - xa)' * rotation_matrix(qa) * Diagonal([1.0 / ax[1]^2, 1.0 / ay[1]^2, 1.0 / az[1]^2]) * rotation_matrix(qa)' * (pa - xa));
        sb[1] - (1.0^2 - (pb - xb)' * rotation_matrix(qb) * Diagonal([1.0 / bx[1]^2, 1.0 / by[1]^2, 1.0 / bz[1]^2]) * rotation_matrix(qb)' * (pb - xb));
    ]
end

function lagrangian(z, y, θ)
    # initialize
    L = 0.0

    # primals
    pa, pb, sa, sb = unpack_primals(z)

    # duals
    λa, λb, νa, νb = unpack_duals(y)

    # objective
    J = objective(z)
    L += J

    # constraints
    c = constraints(z, θ)
    L += sum([λa[1] λb[1]] * c)

    # inequalities
    L -= sa[1] * νa[1]
    L -= sb[1] * νb[1]

    return L
end

L = lagrangian(z, y, θ)
Lz = Symbolics.gradient(L, z)
function residual(w, θ, κ)
    # primals
    z = w[1:8]
    pa, pb, sa, sb = unpack_primals(z)

    # duals
    y = w[8 .+ (1:4)]
    λa, λb, νa, νb = unpack_duals(y)

    # Lagrangian
    lag = lagrangian(z, y, θ)
    lagz = Symbolics.gradient(lag, z)
    con = constraints(z, θ)
    comp = [
            sa[1] * νa[1];
            sb[1] * νb[1];
           ]

    res = [
            lagz;
            con;
            comp .- κ;
          ]

    return res
end

w = [z; y]
r = residual([z; y], θ, κ)
rw = Symbolics.jacobian(r, [z; y])
rθ = Symbolics.jacobian(r, θ)

r_func = eval(Symbolics.build_function(r, w, θ, κ)[2])
rw_func = eval(Symbolics.build_function(rw, w, θ)[2])
rθ_func = eval(Symbolics.build_function(rθ, w, θ)[2])

# pre-allocate
r0 = zeros(nw)
rw0 = zeros(nw, nw)
rθ0 = zeros(nw, nθ)

# random
w0 = randn(nw)
θ0 = randn(nθ)
κ0 = randn(1)

r_func(r0, w0, θ0, κ0)
rw_func(rw0, w0, θ0)
rθ_func(rθ0, w0, θ0)

function rw_func_reg(rw, w, θ)
    rw_func(rw, w, θ)
    rw .+= Diagonal([1.0e-5 * ones(8); -1.0e-5 * ones(4)])
    return
end

(rw0 + Diagonal([1.0e-5 * ones(8); -1.0e-5 * ones(4)])) \ rθ0

## setup
xa = [0.0; 0.0; 0.0]
qa = Quaternion(RotY(0.0) * RotX(0.0 * π))
qa = [qa.w; qa.x; qa.y; qa.z]
ax = 0.1
ay = 0.5
az = 0.1

xb = [1.0; 0.0; 0.0]
qb = Quaternion(RotZ(0.25 * π) * RotY(0.2 * π) * RotX(0.0))
qb = [qb.w; qb.x; qb.y; qb.z]
bx = 0.35
by = 0.1
bz = 0.1

## initialization
w0 = [xa; xb; 1.0; 1.0; 0.0; 0.0; 1.0; 1.0]
θ0 = [xa; qa; xb; qb; ax; ay; az; bx; by; bz]

r_func(r0, w0, θ0, κ0)
rw_func_reg(rw0, w0, θ0)
# rw_func(rw0, w0, θ0)
rθ_func(rθ0, w0, θ0)

rw0 \ rθ0

## solver
idx = RoboDojo.IndicesOptimization(
    nw, nw,
    [[7, 8], [11, 12]], [[7, 8], [11, 12]],
    Vector{Vector{Vector{Int}}}(), Vector{Vector{Vector{Int}}}(),
    collect(1:10),
    collect(11:12),
    Vector{Int}(),
    Vector{Vector{Int}}(),
    collect(11:12),
)

ip = RoboDojo.interior_point(w0, θ0;
    s = RoboDojo.Euclidean(length(w0)),
    idx = idx,
    r! = r_func, rz! = rw_func_reg, rθ! = rθ_func,
    r  = zeros(idx.nΔ),
    rz = zeros(idx.nΔ, idx.nΔ),
    rθ = zeros(idx.nΔ, length(θ0)),
    opts = RoboDojo.InteriorPointOptions(
            undercut=Inf,
            γ_reg=0.1,
            r_tol=1e-6,
            κ_tol=1e-6,
            max_ls=25,
            ϵ_min=0.25,
            diff_sol=true,
            verbose=true))

RoboDojo.interior_point_solve!(ip)

pa = ip.z[0 .+ (1:3)]
pb = ip.z[3 .+ (1:3)]

∂pa∂xa = ip.δz[0 .+ (1:3), 0  .+ (1:3)]
∂pa∂qa = ip.δz[0 .+ (1:3), 3  .+ (1:4)]
∂pa∂xb = ip.δz[0 .+ (1:3), 7  .+ (1:3)]
∂pa∂qb = ip.δz[0 .+ (1:3), 10 .+ (1:4)]

∂pb∂xa = ip.δz[3 .+ (1:3), 0  .+ (1:3)]
∂pb∂qa = ip.δz[3 .+ (1:3), 3  .+ (1:4)]
∂pb∂xb = ip.δz[3 .+ (1:3), 7  .+ (1:3)]
∂pb∂qb = ip.δz[3 .+ (1:3), 10 .+ (1:4)]

# sphere a
sa = GeometryBasics.Sphere(Point(0.0, 0.0, 0.0), 1.0)
color1 = Colors.RGBA(0.7, 0.7, 0.7, 0.5);
setobject!(vis[:spherea], sa, MeshPhongMaterial(color=color1))
settransform!(vis[:spherea], compose(Translation(xa...), LinearMap(Quaternion(qa...) * Diagonal([ax, ay, az]))))

# sphere b
sb = GeometryBasics.Sphere(Point(0.0, 0.0, 0.0), 1.0)
color2 = Colors.RGBA(0.7, 0.7, 0.7, 0.5);
setobject!(vis[:sb], sb, MeshPhongMaterial(color=color2))
settransform!(vis[:sb], compose(Translation(xb...), LinearMap(Quaternion(qb...) * Diagonal([bx, by, bz]))))

# closest points
color_cp = Colors.RGBA(0.0, 0.0, 0.0, 1.0);
cs1 = GeometryBasics.Sphere(Point(0.0, 0.0, 0.0), 0.025)
cs2 = GeometryBasics.Sphere(Point(0.0, 0.0, 0.0), 0.025)
setobject!(vis[:cs1], cs1, MeshPhongMaterial(color=color_cp))
setobject!(vis[:cs2], cs2, MeshPhongMaterial(color=color_cp))
settransform!(vis[:cs1], Translation(pa...))
settransform!(vis[:cs2], Translation(pb...))

ext: jl


text:
function fig_print(x_data, y_data, fig_filename,
					legend_string = nothing, xlabel_string = nothing, ylabel_string = nothing;
					fig_handle = nothing,
					lw_val = 1.5, N_markers = 10, mks_val = 4,
					gfs_val = 9, lfs_val = 8, ar_val = :auto, save_file = 1)
"""
Figure generation and exporting to file

fig_print(x_data, y_data, fig_filename, legend_string, xlabel_string, ylabel_string)

To just save an existing plot "f" to file without modification, use the following code:

fig_print([], [], fig_filename; fig_handle = f)

-------------------
x_data: 1D/2D array, columns are series
y_data: 2D array, 	 columns are series
legend_string: row vector (1D array)

if     fig_handle is given  , then fig_print updates an existing plot
elseif fig_handle is nothing, then fig_print generates a new plot

if 	   x_data or y_data is empty, then fig_print still controls axes, legends, file saving
"""

#------ linespec
linestyle_preset   = [:solid	:solid	:solid	:solid	:dashdot :dash	:dot]
# choose from [:auto, :solid, :dash, :dot, :dashdot, :dashdotdot]
markershape_preset = [:circle	:rect	:cross	:none	:none	 :none	:none]
# choose from [:none, :auto, :circle, :rect, :star5, :diamond, :hexagon, :cross, :xcross, :utriangle, :dtriangle, :rtriangle, :ltriangle, :pentagon, :heptagon, :octagon, :star4, :star6, :star7, :star8, :vline, :hline, :+, :x].

if ~isempty(x_data) && ~isempty(y_data)
#------ marker index
	marker_idx = zeros(size(x_data, 1))
	marker_idx[round.(Int, LinRange(1, size(x_data, 1), N_markers))] = mks_val*ones(N_markers)
	# marker_idx = zeros(length(fig_handle.series_list[1].plotattributes[:x]))
	# marker_idx[round.(Int, LinRange(1,length(fig_handle.series_list[1].plotattributes[:x]), N_markers))] = mks_val*ones(N_markers)

#------ figure generation / augmentation
	if isnothing(fig_handle)
	fig_handle = plot(x_data, y_data, linewidth = lw_val, linestyle = linestyle_preset,
				label = legend_string, legend = :best,
				markersize = marker_idx, markershape = markershape_preset)
	else
		plot!(fig_handle, x_data, y_data, linewidth = lw_val, linestyle = linestyle_preset[:,length(fig_handle.series_list)+1:end],
				label = legend_string, legend = :best,
				markersize = marker_idx, markershape = markershape_preset[:,length(fig_handle.series_list)+1:end])
	end
end

if ~isnothing(fig_handle)
#------ axes
	plot!(fig_handle, aspect_ratio = ar_val)
	plot!(fig_handle, guidefontsize = gfs_val)
	if ~isnothing(xlabel_string) && ~isempty(xlabel_string)
		plot!(fig_handle, xlabel = xlabel_string)
	end
	if ~isnothing(ylabel_string) && ~isempty(ylabel_string)
		plot!(fig_handle, ylabel = ylabel_string)
	end

#------ legend
	if ~isnothing(legend_string) && length(legend_string) > 1 && typeof(legend_string) != String
		for i in 1:length(fig_handle.series_list)
			fig_handle.series_list[i].plotattributes[:label] = legend_string[i]
		end
	end
	plot!(fig_handle, legendfontsize = lfs_val)

#------ save
	if save_file == 1
		fig_dir = "Figures"
		mkpath(fig_dir)
		savefig(fig_handle, joinpath(fig_dir, string("Fig_", fig_filename, ".pdf")))
	end
end

return fig_handle
end


text:
using BenchmarkTools
using Revise
using Base58
using Base64

function int2bytes(x::Integer)
    hex = string(x, base=16)
    if mod(length(hex), 2) != 0
        hex = string("0", hex)
    end
    return hex2bytes(hex)
end

for e ∈ [16,32,64,128,256,512,1024,2048,4096]

    SAMPLE_SIZE = 32

    encode_data = []
    while length(encode_data) < SAMPLE_SIZE
        push!(encode_data, int2bytes(rand(big.(1:big(2)^e))))
    end

    decode_data = []

    i = 1
    encode_avg = 0
    while i < length(encode_data)
        original = @elapsed base58encode(encode_data[i])
        alternative = @elapsed Base58.alt.base58encode(encode_data[i])
        push!(decode_data, alt.base58encode(encode_data[i]))
        encode_avg += original - alternative
        i += 1
    end
    println(length(encode_data[i]) * 8, " bits average encoding Δ: ", round(encode_avg / length(encode_data) * 10^6), " ns")

    i = 1
    decode_avg = 0
    while i < length(decode_data)
        original = @elapsed base58decode(decode_data[i])
        alternative = @elapsed Base58.alt.base58decode(decode_data[i])
        decode_avg += original - alternative
        i += 1
    end
    println(length(encode_data[i]) * 8, " bits average decoding Δ: ", round(decode_avg / length(decode_data) * 10^6), " ns")

end


text:
import Base: show

export YES, NO, nil
export hostname
export release, autorelease # TODO: remove
export loadbundle, framework

@classes NSString, NSNumber, NSDate, NSBundle, NSArray, NSObject, NSHost
export   NSString, NSNumber, NSDate, NSBundle, NSArray, NSObject, NSHost

const YES = true
const NO  = false
const nil = C_NULL

toobject(s::String) = @objc [[NSString alloc] initWithUTF8String:s]

# Conversion of numbers to NSNumber
# See: https://developer.apple.com/documentation/foundation/nsnumber?language=objc
toobject(num::Signed) = @objc [NSNumber numberWithInt:num]
toobject(num::Unsigned) = @objc [NSNumber numberWithUnsignedInt:num]
toobject(num::Float32) = @objc [NSNumber numberWithFloat:num]
toobject(num::Float64) = @objc [NSNumber numberWithDouble:num]
toobject(num::Bool) = @objc [NSNumber numberWithBool:num]

function hostname()
  cstr = @objc [[[NSHost currentHost] localizedName] UTF8String]
  unsafe_string(cstr)
end

release(obj) = @objc [obj release]
autorelease(obj) = @objc [obj autorelease]

# function Base.gc(obj::Object)
#   finalizer(obj, release)
#   obj
# end

function loadbundle(path)
  bundle = @objc [NSBundle bundleWithPath:path]
  bundle.ptr |> Int |> int2bool || error("Bundle $path not found")
  loadedStuff = @objc [bundle load]
  loadedStuff |> int2bool || error("Couldn't load bundle $path")
  return
end

framework(name) = loadbundle("/System/Library/Frameworks/$name.framework")

function show(io::IO, obj::Object)
    addr = repr(UInt(obj.ptr))
    println(io, class(obj), " Object $addr")

    description = @objc [obj description]
    description == nil && return

    description = @objc [description UTF8String]
    description == nil && return

    print(io, unsafe_string(description))
end


text:
using Dates, DataFrames, HDF5, JLD, Requests

include(Pkg.dir("Quandl/src/dataframe.jl"))

dfa = dataframe(load(Pkg.dir("Quandl/test/response.jld"))["asc"])
dfd = dataframe(load(Pkg.dir("Quandl/test/response.jld"))["des"])

facts("dataframe works on Request object") do

  context("there are three rows") do
      @fact  size(dfa,1) => 3
      @fact  size(dfd,1) => 3
  end

  context("there are six columns") do
      @fact  size(dfa,2) => 6
      @fact  size(dfd,2) => 6
  end

  context("correct value at first row, first column, depending on ordering") do
      @fact  dfa[1,2] => 103.02
      @fact  dfd[1,2] => 104.08
  end

  context("oldest date last when ascending") do
      @fact dfa[1,1] < dfa[3,1] => true
  end

  context("oldest date first when descending") do
      @fact dfd[1,1] > dfd[3,1] => true
  end
end


text:
# private helper for union
function _union(A::Multiset{S}, B::Multiset{S}) where {S}
    M = Multiset{S}()
    for (x, v) in A.data
        M[x] = max(v, B[x])
    end
    for (y, w) in B.data
        M[y] = max(A[y], w)
    end
    return M
end

function type_convert(T, A::Multiset)
    M = Multiset{T}()
    for (x, v) in A.data
        xx = convert(T, x)
        M[xx] = v
    end
    return M
end

"""
`union(A,B)` for multisets creates a new multiset in which the
multiplicity of `x` is `max(A[x],B[x])`. This may be invoked as
`A ∪ B`.
"""
function union(A::Multiset{S}, B::Multiset{T}) where {S,T}
    if S == T
        return _union(A, B)
    end
    ST = typejoin(S, T)
    AA = type_convert(ST, A)
    BB = type_convert(ST, B)
    return _union(AA, BB)
end

"""
`A+B` for multisets is the disjoint union, i.e., a new multiset in which the
multiplicity of `x` is `A[x]+B[x]`.
"""
function (+)(A::Multiset{S}, B::Multiset{T}) where {S,T}
    ST = typejoin(S, T)
    M = Multiset{ST}()
    for (x, v) in A.data
        push!(M, x, v)
    end
    for (x, v) in B.data
        push!(M, x, v)
    end
    return M
end

# private helper for A-B
function _minus(A::Multiset{S}, B::Multiset{S}) where {S}
    M = Multiset{S}()
    for (x, v) in A.data
        M.data[x] = max(v - B[x], 0)
    end
    return M
end

"""
`A-B` for multisets is the multiset difference, i.e., a new multiset
in which the multiplicity of `x` is `A[x]-B[x]` unless this goes
below `0`, in which case the multiplicity is 0.
"""
function (-)(A::Multiset{S}, B::Multiset{T}) where {S,T}
    if S == T
        return _minus(A, B)
    end
    ST = typejoin(S, T)
    AA = type_convert(ST, A)
    BB = type_convert(ST, B)
    return _minus(AA, BB)
end

# private helper for intersect
function _intersect(A::Multiset{S}, B::Multiset{S}) where {S}
    M = Multiset{S}()
    for (x, v) in A.data
        push!(M, x, min(v, B[x]))
    end
    return M
end

"""
`intersect(A,B)` for multisets creates a new multiset in which the
multiplicity of `x` is `min(A[x],B[x])`.
This may be abbreviated `A ∩ B`.
"""
function intersect(A::Multiset{S}, B::Multiset{T}) where {S,T}
    if S == T
        return _intersect(A, B)
    end
    ST = typejoin(S, T)
    AA = type_convert(ST, A)
    BB = type_convert(ST, B)
    return _intersect(AA, BB)
end

"""
`A*B` for the Cartesian product of multisets `A` and `B`.
"""
function (*)(A::Multiset{S}, B::Multiset{T}) where {S,T}
    ST = Tuple{S,T}
    M = Multiset{ST}()
    for (a, v) in A.data
        for (b, w) in B.data
            M.data[(a, b)] = v * w
        end
    end
    return M
end

"""
`n*A` is the scalar multiple of a multiset in which the multiplicity of
`x` is `n*A[x]`. Of course, we require `n >= 0`.
"""
function (*)(n::Int, A::Multiset{T}) where {T}
    @assert n >= 0 "Scalar multiplication of a multiset must be by a nonnegative integer"
    M = Multiset{T}()
    for (x, v) in A.data
        M.data[x] = n * v
    end
    return M
end

# private helper for issubset
function _sub(A::Multiset{S}, B::Multiset{S}) where {S}
    for (x, v) in A.data
        if v > B[x]
            return false
        end
    end
    return true
end

function issubset(A::Multiset{S}, B::Multiset{T}) where {S,T}
    if S == T
        return _sub(A, B)
    end
    ST = typejoin(S, T)
    AA = type_convert(ST, A)
    BB = type_convert(ST, B)
    return _sub(AA, BB)
end

(==)(A::Multiset, B::Multiset) = (length(A) == length(B)) && issubset(A, B)

## mixing Multisets with other AbstractSets and the like

union(A::Multiset, B::T) where {T} = union(A, Multiset(B))
union(B::Set, A::Multiset) = union(A, Multiset(B))
union(B::BitSet, A::Multiset) = union(A, Multiset(B))

intersect(A::Multiset, B::T) where {T} = intersect(A, Multiset(B))
intersect(B::Set, A::Multiset) = intersect(A, Multiset(B))
intersect(B::BitSet, A::Multiset) = intersect(A, Multiset(B))

(+)(A::Multiset, B::T) where {T} = A + Multiset(B)
(+)(B::T, A::Multiset) where {T} = A + Multiset(B)

(-)(A::Multiset, B::T) where {T} = A - Multiset(B)
(-)(B::T, A::Multiset) where {T} = Multiset(B) - A

(*)(A::Multiset, B::T) where {T} = A * Multiset(B)
(*)(A::T, B::Multiset) where {T} = Multiset(A) * B


text:
import Pkg
cd("../")
include(joinpath(pwd(), "julenvnew.jl"))
Pkg.activate("GTheoryJulEnvNew")
sep = "/"
working_path = pwd()
settings_path = joinpath(working_path, "settings.yml")
testpath = joinpath(pwd(), "src")
push!(LOAD_PATH, testpath)
using TestOptJulia
using YAML
settings = YAML.load(open(settings_path))
output_path = joinpath(working_path, "output")
start_of_opt(settings, sep, output_path)

max_stars_repo_name: sambuddhac/TestOptJulia

max_issues_repo_name: sambuddhac/TestOptJulia

max_forks_repo_name: sambuddhac/TestOptJulia


text:
## Environment and packages
cd(@__DIR__)
using Pkg; Pkg.activate("."); Pkg.instantiate()

using OrdinaryDiffEq
using ModelingToolkit
using DataDrivenDiffEq
using LinearAlgebra, Optim
using DiffEqFlux, Flux
using Plots
gr()
using JLD2, FileIO
using Statistics
using DelimitedFiles
# Set a random seed for reproduceable behaviour
using Random
Random.seed!(5443)

#### NOTE
# Since the recent release of DataDrivenDiffEq v0.6.0 where a complete overhaul of the optimizers took
# place, SR3 has been used. Right now, STLSQ performs better and has been changed.
# Additionally, the behaviour of the optimization has changed slightly. This has been adjusted
# by decreasing the tolerance of the gradient.

svname = "HudsonBay"
## Data Preprocessing
# The data has been taken from https://jmahaffy.sdsu.edu/courses/f00/math122/labs/labj/q3v1.htm
# Originally published in E. P. Odum (1953), Fundamentals of Ecology, Philadelphia, W. B. Saunders
hudson_bay_data = readdlm("hudson_bay_data.dat", '\t', Float32, '\n')
# Measurements of prey and predator
Xₙ = Matrix(transpose(hudson_bay_data[:, 2:3]))
t = hudson_bay_data[:, 1] .- hudson_bay_data[1, 1]
# Normalize the data; since the data domain is strictly positive
# we just need to divide by the maximum
xscale = maximum(Xₙ, dims =2)
Xₙ .= 1f0 ./ xscale .* Xₙ
# Time from 0 -> n
tspan = (t[1], t[end])

# Plot the data
scatter(t, transpose(Xₙ), xlabel = "t", ylabel = "x(t), y(t)")
plot!(t, transpose(Xₙ), xlabel = "t", ylabel = "x(t), y(t)")

## Direct Identification via SINDy + Collocation

# Create the problem using a gaussian kernel for collocation
full_problem = ContinuousDataDrivenProblem(Xₙ, t, DataDrivenDiffEq.GaussianKernel())
# Look at the collocation
plot(full_problem.t, full_problem.X')
plot(full_problem.t, full_problem.DX')

# Create a Basis
@variables u[1:2]

# Generate the basis functions, multivariate polynomials up to deg 5
# and sine
b = [polynomial_basis(u, 5); sin.(u)]
basis = Basis(b, u)

# Create the thresholds which should be used in the search process
λ = Float32.(exp10.(-7:0.1:5))
# Create an optimizer for the SINDy problem
opt = STLSQ(λ)

# Best result so far
full_res = solve(full_problem, basis, opt, maxiter = 10000, progress = true, denoise = true, normalize = true)

println(full_res)
println(result(full_res))

## Define the network
# Gaussian RBF as activation
rbf(x) = exp.(-(x.^2))

# Define the network 2->5->5->5->2
U = FastChain(
    FastDense(2,5,rbf), FastDense(5,5, rbf), FastDense(5,5, tanh), FastDense(5,2)
)

# Get the initial parameters, first two is linear birth / decay of prey and predator
p = [rand(Float32,2); initial_params(U)]

# Define the hybrid model
function ude_dynamics!(du,u, p, t)
    û = U(u, p[3:end]) # Network prediction
    # We assume a linear birth rate for the prey
    du[1] = p[1]*u[1] + û[1]
    # We assume a linear decay rate for the predator
    du[2] = -p[2]*u[2] + û[2]
end

# Define the problem
prob_nn = ODEProblem(ude_dynamics!,Xₙ[:, 1], tspan, p)

## Function to train the network
# Define a predictor
function predict(θ, X = Xₙ[:,1], T = t)
    Array(solve(prob_nn, Vern7(), u0 = X, p=θ,
                tspan = (T[1], T[end]), saveat = T,
                abstol=1e-6, reltol=1e-6,
                sensealg = ForwardDiffSensitivity()
                ))
end

# Define parameters for Multiple Shooting
group_size = 5
continuity_term = 200.0f0

function loss(data, pred)
	return sum(abs2, data - pred)
end

function shooting_loss(p)
    return multiple_shoot(p, Xₙ, t, prob_nn, loss, Vern7(),
                          group_size; continuity_term)
end

function loss(θ)
    X̂ = predict(θ)
    sum(abs2, Xₙ - X̂) / size(Xₙ, 2) + convert(eltype(θ), 1e-3)*sum(abs2, θ[3:end]) ./ length(θ[3:end])
end

# Container to track the losses
losses = Float32[]

# Callback to show the loss during training
callback(θ,args...) = begin
	l = loss(θ) # Equivalent L2 loss
    push!(losses, l)
    if length(losses)%5==0
        println("Current loss after $(length(losses)) iterations: $(losses[end])")
    end
    false
end

## Training -> First shooting / batching to get a rough estimate

# First train with ADAM for better convergence -> move the parameters into a
# favourable starting positing for BFGS
res1 = DiffEqFlux.sciml_train(shooting_loss, p, ADAM(0.1f0), cb=callback, maxiters = 100)
println("Training loss after $(length(losses)) iterations: $(losses[end])")
# Train with BFGS to achieve partial fit of the data
res2 = DiffEqFlux.sciml_train(shooting_loss, res1.minimizer, BFGS(initial_stepnorm=0.01f0), cb=callback, maxiters = 500)
println("Training loss after $(length(losses)) iterations: $(losses[end])")
# Full L2-Loss for full prediction
res3 = DiffEqFlux.sciml_train(loss, res2.minimizer, BFGS(initial_stepnorm=0.01f0), cb=callback, maxiters = 10000)
println("Final training loss after $(length(losses)) iterations: $(losses[end])")

pl_losses = plot(1:101, losses[1:101], yaxis = :log10, xaxis = :log10, xlabel = "Iterations", ylabel = "Loss", label = "ADAM (Shooting)", color = :blue)
plot!(102:302, losses[102:302], yaxis = :log10, xaxis = :log10, xlabel = "Iterations", ylabel = "Loss", label = "BFGS (Shooting)", color = :red)
plot!(302:length(losses), losses[302:end], color = :black, label = "BFGS (L2)")
savefig(pl_losses, joinpath(pwd(), "plots", "$(svname)_losses.pdf"))

# Rename the best candidate
p_trained = res3.minimizer

## Analysis of the trained network
# Interpolate the solution
tsample = t[1]:0.5:t[end]
X̂ = predict(p_trained, Xₙ[:,1], tsample)
# Trained on noisy data vs real solution
pl_trajectory = scatter(t, transpose(Xₙ), color = :black, label = ["Measurements" nothing], xlabel = "t", ylabel = "x(t), y(t)")
plot!(tsample, transpose(X̂), color = :red, label = ["UDE Approximation" nothing])
savefig(pl_trajectory, joinpath(pwd(), "plots", "$(svname)_trajectory_reconstruction.pdf"))

# Neural network guess
Ŷ = U(X̂,p_trained[3:end])

pl_reconstruction = scatter(tsample, transpose(Ŷ), xlabel = "t", ylabel ="U(x,y)", color = :red, label = ["UDE Approximation" nothing])
plot!(tsample, transpose(Ŷ), color = :red, lw = 2, style = :dash, label = [nothing nothing])
savefig(pl_reconstruction, joinpath(pwd(), "plots", "$(svname)_missingterm_reconstruction.pdf"))
pl_missing = plot(pl_trajectory, pl_reconstruction, layout = (2,1))
savefig(pl_missing, joinpath(pwd(), "plots", "$(svname)_reconstruction.pdf"))
## Symbolic regression via sparse regression (SINDy based)
# We reuse the basis and optimizer defined at the beginning

nn_problem = ContinuousDataDrivenProblem(X̂, tsample, DX = Ŷ)
nn_res = solve(nn_problem, basis, opt, maxiter = 10000, progress = true, normalize = false, denoise = true)
println(nn_res)
println(result(nn_res))

# Define the recovered, hyrid model with the rescaled dynamics
function recovered_dynamics!(du,u, p, t)
    û = nn_res(u, p[3:end]) # Network prediction
    du[1] = p[1]*u[1] + û[1]
    du[2] = -p[2]*u[2] + û[2]
end

p_model = [p_trained[1:2];parameters(nn_res)]

estimation_prob = ODEProblem(recovered_dynamics!, Xₙ[:, 1], tspan, p_model)
# Convert for reuse
sys = modelingtoolkitize(estimation_prob);
dudt = ODEFunction(sys);
estimation_prob = ODEProblem(dudt,Xₙ[:, 1], tspan, p_model)
estimate = solve(estimation_prob, Tsit5(), saveat = t)

##  Fit the found model
function loss_fit(θ)
    X̂ = Array(solve(estimation_prob, Tsit5(), p = θ, saveat = t))
    sum(abs2, X̂ .- Xₙ)
end

# Post-fit the model
res_fit = DiffEqFlux.sciml_train(loss_fit, p_model, BFGS(initial_stepnorm = 0.1f0), maxiters = 1000)
p_fitted = res_fit.minimizer

# Estimate
estimate_rough = solve(estimation_prob, Tsit5(), saveat = 0.1*mean(diff(t)), p = p_model)
estimate = solve(estimation_prob, Tsit5(), saveat = 0.1*mean(diff(t)), p = p_fitted)

# Plot
pl_fitted = plot(t, transpose(Xₙ), style = :dash, lw = 2,color = :black, label = ["Measurements" nothing], xlabel = "t", ylabel = "x(t), y(t)")
plot!(estimate_rough, color = :red, label = ["Recovered" nothing])
plot!(estimate, color = :blue, label = ["Recovered + Fitted" nothing])
savefig(pl_fitted,joinpath(pwd(),"plots","$(svname)recovery_fitting.pdf"))

## Simulation

# Look at long term prediction
t_long = (0.0f0, 50.0f0)
estimate_long = solve(estimation_prob, Tsit5(), saveat = 0.25f0, tspan = t_long,p = p_fitted)
plot(estimate_long.t, transpose(xscale .* estimate_long[:,:]), xlabel = "t", ylabel = "x(t),y(t)")

## Save the results
save(joinpath(pwd(),"results","Hudson_Bay_recovery.jld2"),
    "X", Xₙ, "t" , t, "neural_network" , U, "initial_parameters", p, "trained_parameters" , p_trained, # Training
    "losses", losses, "result", nn_res, "recovered_parameters", parameters(nn_res), # Recovery
    "model", recovered_dynamics!, "model_parameter", p_model, "fitted_parameter", p_fitted,
    "long_estimate", estimate_long) # Estimation

## Post Processing and Plots

c1 = 3 # RGBA(174/255,192/255,201/255,1) # Maroon
c2 = :orange # RGBA(132/255,159/255,173/255,1) # Red
c3 = :blue # RGBA(255/255,90/255,0,1) # Orange
c4 = :purple # RGBA(153/255,50/255,204/255,1) # Purple

p3 = scatter(t, transpose(Xₙ), color = [c1 c2], label = ["x data" "y data"],
             title = "Recovered Model from Hudson Bay Data",
             titlefont = "Helvetica", legendfont = "Helvetica",
             markersize = 5)

plot!(p3,estimate_long, color = [c3 c4], lw=1, label = ["Estimated x(t)" "Estimated y(t)"])
plot!(p3,[19.99,20.01],[0.0,maximum(Xₙ)*1.25],lw=1,color=:black, label = nothing)
annotate!([(10.0,maximum(Xₙ)*1.25,text("Training \nData",12 , :center, :top, :black, "Helvetica"))])
savefig(p3,joinpath(pwd(),"plots","$(svname)full_plot.pdf"))


text:
module bundle-adjustment-primer

# Write your package code here.

end

max_stars_repo_name: ajahraus/bundle-adjustment-primer

max_issues_repo_name: ajahraus/bundle-adjustment-primer

max_forks_repo_name: ajahraus/bundle-adjustment-primer


text:
## First we treat weighted sampling

export sampling_weights
@ap_interface sampling_weights
@ap_sampling_property sampling_weights

sampling_weights(sstyle::UnweightedDiscreteStyle, ap; options...) =
    BasisFunctions.uniformweights(sampling_grid(sstyle, ap; options...))

sampling_weights(sstyle::NormalizedSampling, ap; options...) =
    BasisFunctions.normalizing_weights(sampling_grid(sstyle, ap; options...), measure(ap))

sampling_weights(sstyle::WeightedSampling, ap; options...) =
    sampling_weights(sstyle, platform(ap), platformparameter(ap), samplingparameter(ap); options...)

## Now we can make a sampling operator

export sampling_operator, dual_sampling_operator
@ap_interface sampling_operator
@ap_sampling_property sampling_operator

@ap_interface dual_sampling_operator
@ap_sampling_property dual_sampling_operator

sampling_operator(ss::UnweightedDiscreteStyle, ap::ApproximationProblem; options...) =
    GridSampling(sampling_grid(ss, ap; options...), coefficienttype(ap))

@deprecate samplingoperator(args...; options...) sampling_operator(args...; options...)
@deprecate dualsamplingoperator(args...; options...) dual_sampling_operator(args...; options...)

function sampling_operator(ss::WeightedDiscreteStyle, ap::ApproximationProblem; options...)
    grid = sampling_grid(ss, ap; options...)
    weights = sampling_weights(ss, ap; options...)
    T = coefficienttype(ap)
    S = GridSampling(grid, T)
    D = DiagonalOperator(dest(S), dest(S), weights)
    D*S
end

sampling_operator(ss::DiscreteGramStyle, ap::ApproximationProblem; options...) =
    ProjectionSampling(dictionary(ap), discretemeasure(ap))

sampling_operator(ss::GenericSamplingStyle, ap::PlatformApproximation; options...) =
    sampling_operator(ss, platform(ap), platformparameter(ap), samplingparameter(ap); options...)

sampling_operator(ss::GramStyle, ap::ApproximationProblem; options...) =
    ProjectionSampling(dictionary(ap), measure(ap))

sampling_operator(samplingstyle::RectangularGramStyle, ap::ApproximationProblem;
            projectiondict, options...) =
    ProjectionSampling(projectiondict, measure(ap))

# TODO: fix, remove reference to azdual
dual_sampling_operator(samplingstyle::GramStyle, ap::ApproximationProblem; options...) =
    ProjectionSampling(
        haskey(options, :dualdict) ? options[:dualdict] : azdual(samplingstyle, ap; options...),
            measure(ap))

# TODO: fix, remove reference to azdual
dual_sampling_operator(samplingstyle::DiscreteGramStyle, ap::ApproximationProblem; options...) =
    ProjectionSampling(
        haskey(options, :dualdict) ? options[:dualdict] : azdual(samplingstyle, ap; options...),
            measure(ap))

function dual_sampling_operator(ss::SamplingStyle, ap::ApproximationProblem; options...)
    S = sampling_operator(ap; options...)
    todual(S)
end

todual(S::BasisFunctions.GridSampling) = S
todual(S::BasisFunctions.ProjectionSampling) = S
todual(S::BasisFunctions.GenericCompositeOperator) = compose(S.operators[1], map(inv, S.operators[2:end])...)


text:
module FileReader

    # Подключааем библиотеку для работы с YAML
    using YAML
    # Подключаем библиотеку для работы с JSON
    using JSON

    # Подключаем библиотеку для работы со временем
    include("./Chronometer.jl")
    import .Chronometer

    """Функция импорта настроек из JSON файла

    file_path: путь до файла.
    """
    function from_json(file_path::String)
        @info Chronometer.message_with_time("Выполняется функция чтения данных из JSON файла...")

        @debug "Параметры вызова функции:"
        @debug "- путь до файла в файловой системе (полный путь)" file_path

        result = nothing

        @debug "Инициализируем результат с пустым значенем" result

        try
            @debug "Пытаемся открыть файл в файловой системе"
            file = open(file_path)

            @debug "Результат открытия файла" file
            @debug "Пытаемся прочитать данные из файла"
            data = read(file, String)

            @debug "Результат чтения данных из файла" data
            @debug "Пытаемся распарсить строку JSON в словарь"
            result = JSON.parse(data)

            @debug "Результат парсинга JSON строки" result
        catch error
            @error Chronometer.message_with_time("Не удалось прочитать JSON файл. Ошибка: ") error
        end

        return result
    end

    """Функция импорта настроек из YAML файла

    file_path: путь до файла.
    """
    function from_yaml(file_path::String)::Union{Dict{String, Any}, Nothing}
        @info Chronometer.message_with_time("Выполняется функция чтения данных из YAML файла...")

        @debug "Параметры вызова функции:"
        @debug "> путь до файла в файловой системе (полный путь)" file_path

        result::Union{Dict{String, Any}, Nothing} = nothing
        @debug "Инициализируем результат с пустым значенем" result

        # Читаем файл если путь существует
        if file_path !== nothing
            try
                @debug "Пытаемся прочитать содержимое файла по пути" file_path
                result = YAML.load_file(file_path; dicttype=Dict{String,Any})
                @debug "Результат чтения содержимого файла" result
            catch error
                @error Chronometer.message_with_time("Произошла ошибка при чтении файла $file_path") error
            end
        end

        return result
    end
end


text:
using Test
using DataFrames,DBInterface
using DuckDB
#include("../src/DuckDB.jl")

@testset "DB Connection" begin
    con = DuckDB.connect(":memory:")
    @test isa(con, Base.RefValue{Ptr{Nothing}})
    DuckDB.disconnect(con)
end

@testset "Test README" begin
    con = DuckDB.connect(":memory:")
    res = DuckDB.execute(con, "CREATE TABLE integers(date DATE, jcol INTEGER)")
    res = DuckDB.execute(
        con,
        "INSERT INTO integers VALUES ('2021-09-27', 4), ('2021-09-28', 6), ('2021-09-29', 8)",
    )
    res = DuckDB.execute(con, "SELECT * FROM integers")
    @test isa(DuckDB.toDataFrame(res), DataFrame)
    @test isa(DuckDB.toDataFrame(con, "SELECT * FROM integers"), DataFrame)
    DuckDB.disconnect(con)
end

@testset "HUGE Int test" begin
  con = DuckDB.connect(":memory:")
  res = DuckDB.execute(con,"CREATE TABLE huge(id INTEGER,data HUGEINT);")
  res = DuckDB.execute(con,"INSERT INTO huge VALUES (1,NULL), (2, 1761718171), (3, 171661889178);")
  res = DuckDB.toDataFrame(con,"SELECT * FROM huge")
  DuckDB.disconnect(con)
end

@testset "Interval type" begin
    con = DuckDB.connect(":memory:")
    res = DuckDB.execute(con, "CREATE TABLE interval(interval INTERVAL);")
    res = DuckDB.execute(
        con,
        """
INSERT INTO interval VALUES
(INTERVAL 5 HOUR),
(INTERVAL 12 MONTH),
(INTERVAL 12 MICROSECOND),
(INTERVAL 1 YEAR);
""",
    )
    res = DuckDB.toDataFrame(con, "SELECT * FROM interval;")
    @test isa(res, DataFrame)
    DuckDB.disconnect(con)
end

@testset "Interval type" begin
    con = DuckDB.connect(":memory:")
    res = DuckDB.execute(con, "CREATE TABLE timestamp(timestamp TIMESTAMP , data INTEGER);")
    res = DuckDB.execute(
        con,
        "INSERT INTO timestamp VALUES ('2021-09-27 11:30:00.000', 4), ('2021-09-28 12:30:00.000', 6), ('2021-09-29 13:30:00.000', 8);",
    )
    res = DuckDB.execute(con, "SELECT * FROM timestamp;")
    res = DuckDB.toDataFrame(res)
    @test isa(res, DataFrame)
    DuckDB.disconnect(con)
end

@testset "Items table" begin
    con = DuckDB.connect(":memory:")
    res = DuckDB.execute(
        con,
        "CREATE TABLE items(item VARCHAR, value DECIMAL(10,2), count INTEGER);",
    )
    res = DuckDB.execute(
        con,
        "INSERT INTO items VALUES ('jeans', 20.0, 1), ('hammer', 42.2, 2);",
    )
    res = DuckDB.toDataFrame(con, "SELECT * FROM items;")
    @test isa(res, DataFrame)
    DuckDB.disconnect(con)
end

@testset "Integers and dates table" begin
    db = DuckDB.DB()
    res = DBInterface.execute(db, "CREATE TABLE integers(date DATE, data INTEGER);")
    res = DBInterface.execute(
        db,
        "INSERT INTO integers VALUES ('2021-09-27', 4), ('2021-09-28', 6), ('2021-09-29', 8);",
    )
    res = DBInterface.execute(db, "SELECT * FROM integers;")
    res = DuckDB.toDataFrame(res)
    @test isa(res, DataFrame)
    DBInterface.close!(db)
end


text:
export Sca,ScaUnit,ScaNorm,ScaFT,ScaFTEdge, ScaMid
export Ctr,CtrCorner,CtrFFT,CtrFT,CtrMid,CtrEnd, CtrRFFT,CtrRFT  # These do not work, as they need the size information about the real-space array

# define types to specify where the center point is
abstract type Ctr end  # Center of the array
struct CtrCorner <: Ctr end  # corner voxel is zero
struct CtrFFT <: Ctr end # corresponding to FFTs
struct CtrRFFT <: Ctr end # corresponding to RFFTs
struct CtrFT <: Ctr end # corresponding to FTs  (meaning shifted FFTs)
struct CtrRFT <: Ctr end # corresponding to RFTs
struct CtrMid <: Ctr end # middle of the array
struct CtrEnd <: Ctr end # other corner voxel is zero

"""
    Ctr

Abstract type to specify the reference position
from which several other types subtype.

# Possible subtypes
* `CtrCorner`: Set the reference pixel in the corner
* `CtrFFT`: Set the reference pixel to the FFT center.
* `CtrFT`: Set the reference pixel to the FT center. FT means that the zero frequency is at the FFT convention center (`size ÷ 2 + 1`).
* `CtrRFFT`: Set the reference pixel to the RFFT center. Same as `CtrFFT` but the first dimension has center at 1.
* `CtrRFT`: Set the reference pixel to the RFT center. FT means that the zero frequency is at the FFT convention center (`size ÷ 2 + 1`).
            Same as `CtrFT` but the first dimension has center at 1.
* `CtrMid`: Set the reference pixel to real mid. For uneven arrays it is the center pixel, for even arrays it is the centered around a half pixel.
* `CtrEnd` Set the reference to the end corner (last pixel)
"""
Ctr

# In the equations below, it is important to not use 1.0 or 0.0 as these would cause the result to be of type FLoat64
get_offset(size, ::Type{CtrCorner}) = size.*0 .+ 1
get_offset(size, ::Type{CtrFT}) = size.÷2 .+ 1
get_offset(size, ::Type{CtrRFT}) = Base.setindex(size.÷2,0,1) .+ 1
get_offset(size, ::Type{CtrFFT}) = size.*0 .+ 1
get_offset(size, ::Type{CtrRFFT}) = size.*0 .+ 1
get_offset(size, ::Type{CtrMid}) = (size.+1)./2
get_offset(size, ::Type{CtrEnd}) = size.+0
get_offset(size, t::Number) = ntuple(i -> t, length(size))
get_offset(dummy, t::NTuple) = t
get_offset(dummy, t::IterType) = t
get_offset(dummy, t::Matrix) = Tuple(Tuple(t[:,n]) for n in 1:size(t,2))  # converts the matrix to an iterable collection for convenience

abstract type Sca end # scaling of the array
struct ScaUnit <: Sca end # pixel distance is one
struct ScaNorm <: Sca end # total size along each dimension normalized to 1.0
struct ScaMid <: Sca end # such that the border pixel always reaches 1.0
struct ScaFT <: Sca end # reciprocal Fourier coordinates
struct ScaRFT <: Sca end # reciprocal Fourier coordinates for rFTs.
struct ScaFTEdge <: Sca end # such that the edge of the Fourier space is 1.0
struct ScaRFTEdge <: Sca end # such that the edge of the Fourier space is 1.0

"""
    Sca

Abstract type to indicate a scaling from which several other types subtype.

# Possible subtypes
* `ScaUnit`: No scaling of the indices
* `ScaNorm`: Total length along each dimension is normalized to 1
* `ScaMid`: Reaches 1.0 at the borders, if used in combination with `CtrMid`. Useful for keeping real-space symmetry.
* `ScaFT`: Reciprocal Fourier coordinates compared to Nyquist sampling
* `ScaFTEdge`: Such that the edge (in FFT sense) of the pixel is 1.0
"""
Sca

get_scale(sz, ::Type{ScaUnit}) = ntuple(_ -> one(Int), length(sz))
get_scale(sz, ::Type{ScaNorm}) = 1 ./ (max.(sz .- 1, 1))
get_scale(sz, ::Type{ScaMid}) = 1 ./ (max.((sz .- 1)./2, 1)) # reaches 1.0 at the outermost pixel.
get_scale(sz, ::Type{ScaFT}) = 0.5 ./ (max.(sz ./ 2, 1))
# get_scale(size, ::Type{ScaRFT}) = 0.5 ./ (max.(Base.setindex(size./ 2,size[1]-1,1), 1))  # These scales are wrong! They need the information on the real-space size!
get_scale(sz, ::Type{ScaFTEdge}) = 1 ./ (max.(sz ./ 2, 1))
# get_scale(size, ::Type{ScaRFTEdge}) = 1 ./ (max.(Base.setindex(size./ 2,size[1]-1,1), 1))
get_scale(sz, t::Number) = ntuple(i -> t, length(sz))
get_scale(dummy, t::NTuple) = t
get_scale(dummy, t::IterType) = t
get_scale(dummy, t::Matrix) = Tuple(Tuple(t[:,n]) for n in 1:size(t,2))  # converts the matrix to an iterable collection for convenience

"""
    apply_tuple_list(f, t1,t2)
    applies a two-argument function to tubles and iterables of tuples
Example:
```jldoctest
julia> IndexFunArrays.apply_tuple_list.(.*,2.0,((1.0,2.0),(2.0,3.0),(3.0,3.5),(4.0,5.5)))
((2.0, 4.0), (4.0, 6.0), (6.0, 7.0), (8.0, 11.0))
```
"""
function apply_tuple_list(f, t1, t2)  # applies a two-argument function to tubles and iterables of tuples
    return f(t1,t2)
end

function apply_tuple_list(f, t1, t2::Nothing)  # applies a two-argument function to tubles and iterables of tuples
    return t2
end

function apply_tuple_list(f, t1::Nothing, t2)  # applies a two-argument function to tubles and iterables of tuples
    return t1
end

function apply_tuple_list(f, t1, t2::IterType)
    return Tuple([f(t1,a2) for a2 in t2])
end

function apply_tuple_list(f, t1::IterType, t2)
    res= Tuple([f(a1,t2) for a1 in t1])
    return res
end

function apply_tuple_list(f, t1::IterType, t2::IterType)
    return Tuple([f(a[1],a[2]) for a in zip(t1,t2)])
end


text:
using Test
using SAC

@testset "Utility" begin
    # Component flipping
    let s = SACtr(rand(3), 1)
        # Horizontal component pointing at 15°
        s[:cmpaz], s[:cmpinc] = 15, 90
        @test flip_component(s)[:cmpaz] ≈ 195
        @test flip_component(s)[:cmpinc] ≈ 90
        # Arbitrary component
        s[:cmpaz], s[:cmpinc] = 271, 30
        s′ = flip_component!(s)
        @test s′[:cmpaz] ≈ 91
        @test s′[:cmpinc] ≈ 150
        @test flip(s) == flip_component(s)
    end

    # Trace swapping
    let s1 = SACtr(rand(3), 1), s2 = SACtr(rand(3), 1), s1′ = deepcopy(s1), s2′ = deepcopy(s2)
        SAC.swap_traces!(s1, s2)
        @test s1′ == s2
        @test s2′ == s1
    end
end


text:
@testset "Test setup for plotting PhyloNetworks objects" begin

  net = readTopology("(A:2.5,((B:1,#H1:0.5::0.1):1,(C:1,(D:0.5)#H1:0.5::0.9):1):0.5);")
  @test PhyloPlots.getEdgeNodeCoordinates(net, true) == (
      [1., 2.5,2.5,1.5,2.5,3., 2.5,1.5, 1.],
      [3.5,3.5,3., 2.5,3.5,3.5,3., 2.5, 1.5],
      [4., 3., 3., 3., 2., 1., 1., 1.5, 2.25],
      [4., 3., 1., 3., 2., 1., 1., 1.5, 2.25],
      [3.5,3.5,2.5,3.5,3.5,3., 2.5,1.5, 1.],
      [4., 3., 3., 2., 1., 1., 1.5,2.25,3.125],
      [0., 0., 3., 0., 0., 1., 1., 1.5, 2.25],
      [0., 0., 3., 0., 0., 1., 2., 3.,  4.],
      1., 3.5, 1., 4.)
  # in this order:
  # edge_xB, edge_xE, edge_yB, edge_yE,
  # node_x, node_y, node_yB, node_yE,
  # xmin, xmax, ymin, ymax
  @test PhyloPlots.getEdgeNodeCoordinates(net, false) == (
      [1.,3.,3.,2.,3.,4.,3.,2., 1.],
      [5.,5.,4.,3.,5.,5.,4.,3., 2.],
      [4.,3.,3.,3.,2.,1.,1.,1.5,2.25],
      [4.,3.,1.,3.,2.,1.,1.,1.5,2.25],
      [5.,5.,3.,5.,5.,4.,3., 2.,  1.],
      [4.,3.,3.,2.,1.,1.,1.5,2.25,3.125],
      [0.,0.,3.,0.,0.,1.,1., 1.5, 2.25],
      [0.,0.,3.,0.,0.,1.,2., 3.,  4.],
      1., 5.,1., 4.)
  dat = DataFrame(node=[-5,-3,-4,5,100],bs=["90","95","99","mytip","bogus"],edge=[8,9,4,6,200]);
  @test_logs (:warn, "Some node numbers in the nodeLabel data frame are not found in the network:\n 100") PhyloPlots.checkNodeDataFrame(net, dat);
  @test_logs (:warn, "nodeLabel should have 2+ columns, the first one giving the node numbers (Integer)") PhyloPlots.checkNodeDataFrame(net, dat[!,2:3])
  dat = DataFrame(node=[-5,-3,missing,5],
                  bs=["90","95","99","mytip"],edge=[8,9,4,6]);
  @test PhyloPlots.checkNodeDataFrame(net, dat) == (true,
    DataFrame(node=[-5,-3,5],bs=["90","95","mytip"],edge=[8,9,6]))
  dat = DataFrame(node=[-5,-3,-4,5],bs=["90","95","99","mytips"]);
  @test PhyloPlots.prepareNodeDataFrame(net,dat,true,true,true,collect(1.:9),collect(10.:18)) ==
    DataFrame(name=["A","B","","C","D","H1","","",""],
    num=["1","2","-4","4","5","3","-5","-3","-2"],
    lab=["","","99","","mytips","","90","95",""],
    lea=[true,true,false,true,true,false,false,false,false],
    x=collect(1.:9), y=collect(10.:18))
  dat = DataFrame(edge=[8,9,4,6,200],bs=["90","95","99","mytips","bogus"]);
  @test_logs (:warn, "Some edge numbers in the edgeLabel data frame are not found in the network:\n 200") PhyloPlots.prepareEdgeDataFrame(
    net,dat,true,collect(1.:9),collect(10.:18),collect(19.:27),collect(28.:36));
  dat = DataFrame(edge=[8,9,4,6],bs=[missing,"95","99","mytips"]);
  @test PhyloPlots.prepareEdgeDataFrame(net,dat,false,collect(1.:9),collect(11.:19),
    collect(21.:29),collect(31.:39)) == (true, DataFrame(
    len=["2.5","1","0.5","1","1","0.5","0.5","1","0.5"],
    gam=["1","1","0.1","1","1","1","0.9","1","1"],
    num=["1","2","3","4","5","6","7","8","9"],
    lab=["","","","99","","mytips","","","95"], # if not missing: "90" second to last (row 8)
    hyb=[false,false,true,false,false,false,true,false,false],
    min=[false,false,true,false,false,false,false,false,false],
    x=collect(6.:14),y=collect(26.:34)))

  # example with level-2 network, non-tree child:
  # one hybrid node ends up as a leaf in the major tree.
  # no major child edge to follow to set coordinates
  net = readTopology("((((B)#H1:::0.2)#H2,((D,C,#H2:::0.8)S1,(#H1,A)S2)S3)S4);")
  @test_logs plot(net, :R, showNodeNumber=true, showGamma=true);
  @test PhyloPlots.getEdgeNodeCoordinates(net, false) == (
    [5.,4.,1.,3.,3.,3.,2.,4.,4.,2.,1.],
    [6.,5.,4.,6.,6.,4.,3.,5.,6.,4.,2.],
    [2.,2.1,2.275,4.,3.,2.1,3.05,2.,1.,1.5,2.275],
    [2.,2. ,2.1,  4.,3.,2.1,3.05,2.,1.,1.5,2.275],
    [6.,5.,4., 6.,6.,3.,  6.,4., 2.,   1.],
    [2.,2.,2.1,4.,3.,3.05,1.,1.5,2.275,2.275],
    [0.,2.,2.1,0.,0.,2.1,0.,1.,1.5, 2.275],
    [0.,2.,2.1,0.,0.,4. ,0.,2.,3.05,2.275],
    1.,6.,1.,4.)
  net = readTopology("((((B)#H1:::0.2)#H2,((D,C,#H2)S1,(#H1,A)S2)S3)S4);")
  @test_logs plot(net, :R, showNodeNumber=true, showGamma=true);
  @test PhyloPlots.getEdgeNodeCoordinates(net, false) == (
    [5.,4.,1.,3.,3.,3.,2.,4.,4.,2.,1.],
    [6.,5.,4.,6.,6.,4.,3.,5.,6.,4.,2.],
    [2.,2.1,2.1,4.,3.,3.5,3.5,2.,1.,1.5,2.5],
    [2.,2. ,2.1,4.,3.,2.1,3.5,2.,1.,1.5,2.5],
    [6.,5.,4., 6.,6.,3., 6.,4., 2., 1.],
    [2.,2.,2.1,4.,3.,3.5,1.,1.5,2.5,2.3],
    [0.,2.,2.1,0.,0.,3.,0.,1.,1.5,2.1],
    [0.,2.,2.1,0.,0.,4.,0.,2.,3.5,2.5],
    1.,6.,1.,4.0)
end


text:
########################### Interfaces ########################################################
get_variable_key(variabletype, d) = error("Not Implemented")

get_variable_binary(pv, t::Type{<:PSY.Component}, _) =
    error("`get_variable_binary` must be implemented for $pv and $t")

get_variable_warm_start_value(_, ::PSY.Component, __) = nothing

get_variable_lower_bound(_, ::PSY.Component, __) = nothing

get_variable_upper_bound(_, ::PSY.Component, __) = nothing

get_multiplier_value(x, y::PSY.Component, z) =
    error("Unable to get parameter $x for device $y for formulation $z")

get_expression_type_for_reserve(_, y::Type{<:PSY.Component}, z) =
    error("`get_expression_type_for_reserve` must be implemented for $y and $z")

get_initial_conditions_device_model(
    ::OperationModel,
    ::DeviceModel{T, D},
) where {T <: PSY.Device, D <: AbstractDeviceFormulation} =
    error("`get_initial_conditions_device_model` must be implemented for $T and $D")

requires_initialization(::AbstractDeviceFormulation) = false

does_subcomponent_exist(T::PSY.Component, S::Type{<:PSY.Component}) =
    error("`does_subcomponent_exist` must be implemented for $T and subcomponent type $S")

_get_initial_condition_type(
    X::Type{<:ConstraintType},
    Y::Type{<:PSY.Component},
    Z::Type{<:AbstractDeviceFormulation},
) = error("`_get_initial_condition_type` must be implemented for $X , $Y and $Z")

get_initial_conditions_device_model(
    ::OperationModel,
    model::DeviceModel{T, FixedOutput},
) where {T <: PSY.Device} = model


text:
module linearMeshIO

include("./gmsh/GMSHtypes.jl")
include("./openfoam/ofMeshTypes.jl")
include("./tgrid/TGridTypes.jl")
include("./vtk/VTKtypes.jl")

include("./gmsh/GMSHreader.jl")
include("./openfoam/ofMeshReader.jl")
include("./tgrid/TGridReader.jl")
include("./vtk/VTKreader.jl")

end

max_stars_repo_name: c-schubert/linearMeshIO

max_issues_repo_name: c-schubert/linearMeshIO

max_forks_repo_name: c-schubert/linearMeshIO


text:
using Documenter
using Gadfly
using FFAST
using Weave

function addNISTHeaders(htmlfile::String)
    # read HTML
    html = transcode(String,read(htmlfile))
    # Find </head>
    i = findfirst(r"</[Hh][Ee][Aa][Dd]>", html)
    # Already added???
    j = findfirst("nist-header-footer", html)
    if isnothing(j) && (!isnothing(i))
        # Insert nist-pages links right before </head>
        res = html[1:i.start-1]*
            "<link rel=\"stylesheet\" href=\"https://pages.nist.gov/nist-header-footer/css/nist-combined.css\">\n"*
            "<script src=\"https://pages.nist.gov/nist-header-footer/js/jquery-1.9.0.min.js\" type=\"text/javascript\" defer=\"defer\"></script>\n"*
            "<script src=\"https://pages.nist.gov/nist-header-footer/js/nist-header-footer.js\" type=\"text/javascript\" defer=\"defer\"></script>\n"*
            html[i.start:end]
        write(htmlfile, res)
        println("Inserting NIST header/footer into $htmlfile")
    end
    return htmlfile
end

weaveit(name) = weave(joinpath("src", "$name"), out_path=joinpath("src", "$(splitext(name)[1]).md"), doctype="github")

names = ( "example.jmd", )

weaveit.(names)

makedocs(
    modules = [FFAST],
    sitename = "FFAST.jl",
    pages = [ "Home" => "index.md", "Using FFAST" => "example.md" ]
)

map(name->rm(joinpath("src","$(splitext(name)[1]).md")), names)
addNISTHeaders(joinpath(@__DIR__, "build","index.html"))
addNISTHeaders.(map(name->joinpath(@__DIR__, "build", splitext(name)[1], "index.html"), names))

# deploydocs(repo = "github.com/NicholasWMRitchie/FFAST.jl.git")


text:
using MultivariateStats
using Test
using LinearAlgebra
using Random

@testset "Ridge Regression" begin

    Random.seed!(34568)

    ## data

    m = 9
    n = 6
    n2 = 3

    X = randn(m, n)
    A = randn(n, n2)
    Xt = X'

    b = randn(1, n2)

    E = randn(m, n2) * 0.1
    Y0 = X * A + E
    Y1 = X * A .+ b + E

    y0 = Y0[:,1]
    y1 = Y1[:,1]

    ## llsq

    A = llsq(X, Y0; dims=1, bias=false)
    A_r = copy(A)
    @test size(A) == (n, n2)
    @test X'Y0 ≈ X'X * A

    a = llsq(X, y0; dims=1, bias=false)
    @test size(a) == (n,)
    @test a ≈ A[:,1]

    A = llsq(Xt, Y0; dims=2, bias=false)
    @test size(A) == (n, n2)
    @test A ≈ A_r

    a = llsq(Xt, y0; dims=2, bias=false)
    @test size(a) == (n,)
    @test a ≈ A[:,1]

    Aa = llsq(X, Y1; dims=1, bias=true)
    Aa_r = copy(Aa)
    @test size(Aa) == (n+1, n2)
    A, b = Aa[1:end-1,:], Aa[end:end,:]
    @test X' * (Y1 .- b) ≈ X'X * A

    aa = llsq(X, y1; dims=1, bias=true)
    @test aa ≈ Aa[:,1]

    Aa = llsq(Xt, Y1; dims=2, bias=true)
    @test Aa ≈ Aa_r

    aa = llsq(Xt, y1; dims=2, bias=true)
    @test aa ≈ Aa[:,1]

    # test default dims=2 argument
    aa = llsq(X, y1)
    @test aa ≈ Aa[:,1]
    @test_throws DimensionMismatch llsq(X, y1; dims=2)
    @test_throws DimensionMismatch llsq(Xt, y1)

    ## ridge (with Real r)

    r = 0.1

    A = ridge(X, Y0, r; dims=1, bias=false)
    A_r = copy(A)
    @test size(A) == (n, n2)
    @test X'Y0 ≈ (X'X + r * I) * A

    a = ridge(X, y0, r; dims=1, bias=false)
    @test size(a) == (n,)
    @test a ≈ A[:,1]

    A = ridge(Xt, Y0, r; dims=2, bias=false)
    @test size(A) == (n, n2)
    @test A ≈ A_r

    a = ridge(Xt, y0, r; dims=2, bias=false)
    @test size(a) == (n,)
    @test a ≈ A[:,1]

    Aa = ridge(X, Y1, r; dims=1, bias=true)
    Aa_r = copy(Aa)
    @test size(Aa) == (n+1, n2)
    A, b = Aa[1:end-1,:], Aa[end:end,:]
    @test X' * (Y1 .- b) ≈ (X'X + r * I) * A

    aa = ridge(X, y1, r; dims=1, bias=true)
    @test aa ≈ Aa[:,1]

    Aa = ridge(Xt, Y1, r; dims=2, bias=true)
    @test Aa ≈ Aa_r

    aa = ridge(Xt, y1, r; dims=2, bias=true)
    @test aa ≈ Aa[:,1]

    # test default dims=2 argument
    aa = ridge(X, y1, r)
    @test aa ≈ Aa[:,1]
    @test_throws DimensionMismatch ridge(X, y1, r; dims=2)
    @test_throws DimensionMismatch ridge(Xt, y1, r)

    ## ridge (with diagonal r)

    r = 0.05 .+ 0.1 .* rand(n)

    A = ridge(X, Y0, r; dims=1, bias=false)
    A_r = copy(A)
    @test size(A) == (n, n2)
    @test X'Y0 ≈ (X'X + diagm(0=>r)) * A

    a = ridge(X, y0, r; dims=1, bias=false)
    @test size(a) == (n,)
    @test a ≈ A[:,1]

    A = ridge(Xt, Y0, r; dims=2, bias=false)
    @test size(A) == (n, n2)
    @test A ≈ A_r

    a = ridge(Xt, y0, r; dims=2, bias=false)
    @test size(a) == (n,)
    @test a ≈ A[:,1]

    Aa = ridge(X, Y1, r; dims=1, bias=true)
    Aa_r = copy(Aa)
    @test size(Aa) == (n+1, n2)
    A, b = Aa[1:end-1,:], Aa[end:end,:]
    @test X' * (Y1 .- b) ≈ (X'X + diagm(0=>r)) * A

    aa = ridge(X, y1, r; dims=1, bias=true)
    @test aa ≈ Aa[:,1]

    Aa = ridge(Xt, Y1, r; dims=2, bias=true)
    @test Aa ≈ Aa_r

    aa = ridge(Xt, y1, r; dims=2, bias=true)
    @test aa ≈ Aa[:,1]

    ## ridge (with qudratic r matrix)

    Q = qr(randn(n, n)).Q
    r = Q' * diagm(0=>r) * Q

    A = ridge(X, Y0, r; dims=1, bias=false)
    A_r = copy(A)
    @test size(A) == (n, n2)
    @test X'Y0 ≈ (X'X + r) * A

    a = ridge(X, y0, r; dims=1, bias=false)
    @test size(a) == (n,)
    @test a ≈ A[:,1]

    A = ridge(Xt, Y0, r; dims=2, bias=false)
    @test size(A) == (n, n2)
    @test A ≈ A_r

    a = ridge(Xt, y0, r; dims=2, bias=false)
    @test size(a) == (n,)
    @test a ≈ A[:,1]

    Aa = ridge(X, Y1, r; dims=1, bias=true)
    Aa_r = copy(Aa)
    @test size(Aa) == (n+1, n2)
    A, b = Aa[1:end-1,:], Aa[end:end,:]
    @test X' * (Y1 .- b) ≈ (X'X + r) * A

    aa = ridge(X, y1, r; dims=1, bias=true)
    @test aa ≈ Aa[:,1]

    Aa = ridge(Xt, Y1, r; dims=2, bias=true)
    @test Aa ≈ Aa_r

    aa = ridge(Xt, y1, r; dims=2, bias=true)
    @test aa ≈ Aa[:,1]

end


text:
"""
    struct TabularV <: AbstractVApproximator{Int}
        table::Vector{Float64}
    end

Using a `table` of type `Vector{Float64}` to record the state values.
"""
struct TabularV <: AbstractVApproximator{Int}
    table::Vector{Float64}
end

TabularV(ns::Int, init::Float64=0.) = TabularV(fill(init, ns))

(VLearner::TabularV)(s::Int) = VLearner.table[s]

update!(VLearner::TabularV, s::Int, e::Float64) = VLearner.table[s] += e


text:
using FIRRemez
using Test

@testset "FIRRemez.jl" begin
    # Write your tests here.
end

max_stars_repo_name: RoyCCWang/FIRRemez

max_issues_repo_name: RoyCCWang/FIRRemez

max_forks_repo_name: RoyCCWang/FIRRemez


text:
@testset "simple2" begin
    @test true
    push!(Hijack.RUN, 5)
end

max_issues_repo_name: rfourquet/InlineTest

max_forks_repo_name: rfourquet/InlineTest


text:
using Plots
plotlyjs()

include("../julia/pprocess_inhomopoisson.jl");

f = 0.5;
lambda = 20*(1+sin(2.0 * pi * f * ts));
x_inhomopoisson, rs_inhomopoisson = pprocess_inhomopoisson(lambda);

plot(ts, rs_inhomopoisson)
scatter!(x_inhomopoisson, zeros(length(x_inhomopoisson)))
xlabel!("Time (s)")
ylabel!("Firing rate (Hz)")

max_stars_repo_name: ganow/julia-point-process

max_issues_repo_name: ganow/julia-point-process

max_forks_repo_name: ganow/julia-point-process


text:
################################################################################
#                             Projection utilities                             #
################################################################################

function project_position(scene, point, model)

    # use transform func
    point = Makie.apply_transform(scene.transformation.transform_func[], point)

    res = scene.camera.resolution[]
    p4d = to_ndim(Vec4f, to_ndim(Vec3f, point, 0f0), 1f0)
    clip = scene.camera.projectionview[] * model * p4d
    @inbounds begin
    # between -1 and 1
        p = (clip ./ clip[4])[Vec(1, 2)]
        # flip y to match cairo
        p_yflip = Vec2f(p[1], -p[2])
        # normalize to between 0 and 1
        p_0_to_1 = (p_yflip .+ 1f0) / 2f0
    end
    # multiply with scene resolution for final position
    return p_0_to_1 .* res
end

project_scale(scene::Scene, s::Number, model = Mat4f(I)) = project_scale(scene, Vec2f(s), model)

function project_scale(scene::Scene, s, model = Mat4f(I))
    p4d = to_ndim(Vec4f, s, 0f0)
    p = @inbounds (scene.camera.projectionview[] * model * p4d)[Vec(1, 2)] ./ 2f0
    return p .* scene.camera.resolution[]
end

function project_rect(scene, rect::Rect, model)
    mini = project_position(scene, minimum(rect), model)
    maxi = project_position(scene, maximum(rect), model)
    return Rect(mini, maxi .- mini)
end

function project_polygon(scene, poly::P, model) where P <: Polygon
    ext = decompose(Point2f, poly.exterior)
    interiors = decompose.(Point2f, poly.interiors)
    Polygon(
        Point2f.(project_position.(Ref(scene), ext, Ref(model))),
        [Point2f.(project_position.(Ref(scene), interior, Ref(model))) for interior in interiors],
    )
end

scale_matrix(x, y) = Cairo.CairoMatrix(x, 0.0, 0.0, y, 0.0, 0.0)

########################################
#          Rotation handling           #
########################################

function to_2d_rotation(x)
    quat = to_rotation(x)
    return -Makie.quaternion_to_2d_angle(quat)
end

function to_2d_rotation(::Makie.Billboard)
    @warn "This should not be reachable!"
    0
end

remove_billboard(x) = x
remove_billboard(b::Makie.Billboard) = b.rotation

to_2d_rotation(quat::Makie.Quaternion) = -Makie.quaternion_to_2d_angle(quat)

# TODO: this is a hack around a hack.
# Makie encodes the transformation from a 2-vector
# to a quaternion as a rotation around the Y-axis,
# when it should be a rotation around the X-axis.
# Since I don't know how to fix this in GLMakie,
# I've reversed the order of arguments to atan,
# such that our behaviour is consistent with GLMakie's.
to_2d_rotation(vec::Vec2f) = atan(vec[1], vec[2])

to_2d_rotation(n::Real) = n

################################################################################
#                                Color handling                                #
################################################################################

function rgbatuple(c::Colorant)
    rgba = RGBA(c)
    red(rgba), green(rgba), blue(rgba), alpha(rgba)
end

function rgbatuple(c)
    colorant = to_color(c)
    if !(colorant isa Colorant)
        error("Can't convert $(c) to a colorant")
    end
    return rgbatuple(colorant)
end

to_uint32_color(c) = reinterpret(UInt32, convert(ARGB32, c))

function numbers_to_colors(numbers::AbstractArray{<:Number}, primitive)

    colormap = get(primitive, :colormap, nothing) |> to_value |> to_colormap
    colorrange = get(primitive, :colorrange, nothing) |> to_value

    if colorrange === Makie.automatic
        colorrange = extrema(numbers)
    end

    Makie.interpolated_getindex.(
        Ref(colormap),
        Float64.(numbers), # ints don't work in interpolated_getindex
        Ref(colorrange))
end

########################################
#     Image/heatmap -> ARGBSurface     #
########################################

function to_cairo_image(img::AbstractMatrix{<: AbstractFloat}, attributes)
    to_cairo_image(to_rgba_image(img, attributes), attributes)
end

function to_rgba_image(img::AbstractMatrix{<: AbstractFloat}, attributes)
    Makie.@get_attribute attributes (colormap, colorrange, nan_color, lowclip, highclip)

    nan_color = Makie.to_color(nan_color)
    lowclip = isnothing(lowclip) ? lowclip : Makie.to_color(lowclip)
    highclip = isnothing(highclip) ? highclip : Makie.to_color(highclip)

    [get_rgba_pixel(pixel, colormap, colorrange, nan_color, lowclip, highclip) for pixel in img]
end

to_rgba_image(img::AbstractMatrix{<: Colorant}, attributes) = RGBAf.(img)

function get_rgba_pixel(pixel, colormap, colorrange, nan_color, lowclip, highclip)
    vmin, vmax = colorrange

    if isnan(pixel)
        RGBAf(nan_color)
    elseif pixel < vmin && !isnothing(lowclip)
        RGBAf(lowclip)
    elseif pixel > vmax && !isnothing(highclip)
        RGBAf(highclip)
    else
        RGBAf(Makie.interpolated_getindex(colormap, pixel, colorrange))
    end
end

function to_cairo_image(img::AbstractMatrix{<: Colorant}, attributes)
    to_cairo_image(to_uint32_color.(img), attributes)
end

function to_cairo_image(img::Matrix{UInt32}, attributes)
    # we need to convert from column-major to row-major storage,
    # therefore we permute x and y
    return Cairo.CairoARGBSurface(permutedims(img))
end

################################################################################
#                                Mesh handling                                 #
################################################################################

struct FaceIterator{Iteration, T, F, ET} <: AbstractVector{ET}
    data::T
    faces::F
end

function (::Type{FaceIterator{Typ}})(data::T, faces::F) where {Typ, T, F}
    FaceIterator{Typ, T, F}(data, faces)
end
function (::Type{FaceIterator{Typ, T, F}})(data::AbstractVector, faces::F) where {Typ, F, T}
    FaceIterator{Typ, T, F, NTuple{3, eltype(data)}}(data, faces)
end
function (::Type{FaceIterator{Typ, T, F}})(data::T, faces::F) where {Typ, T, F}
    FaceIterator{Typ, T, F, NTuple{3, T}}(data, faces)
end
function FaceIterator(data::AbstractVector, faces)
    if length(data) == length(faces)
        FaceIterator{:PerFace}(data, faces)
    else
        FaceIterator{:PerVert}(data, faces)
    end
end

Base.size(fi::FaceIterator) = size(fi.faces)
Base.getindex(fi::FaceIterator{:PerFace}, i::Integer) = fi.data[i]
Base.getindex(fi::FaceIterator{:PerVert}, i::Integer) = fi.data[fi.faces[i]]
Base.getindex(fi::FaceIterator{:Const}, i::Integer) = ntuple(i-> fi.data, 3)

color_or_nothing(c) = c === nothing ? nothing : to_color(c)

function per_face_colors(
        color, colormap, colorrange, matcap, vertices, faces, normals, uv,
        lowclip=nothing, highclip=nothing, nan_color=nothing
    )
    if matcap !== nothing
        wsize = reverse(size(matcap))
        wh = wsize .- 1
        cvec = map(normals) do n
            muv = 0.5n[Vec(1,2)] .+ Vec2f(0.5)
            x, y = clamp.(round.(Int, Tuple(muv) .* wh) .+ 1, 1, wh)
            return matcap[end - (y - 1), x]
        end
        return FaceIterator(cvec, faces)
    elseif color isa Colorant
        return FaceIterator{:Const}(color, faces)
    elseif color isa AbstractArray
        if color isa AbstractVector{<: Colorant}
            return FaceIterator(color, faces)
        elseif color isa AbstractArray{<: Number}
            low, high = extrema(colorrange)
            cvec = map(color[:]) do c
                if isnan(c) && nan_color !== nothing
                    return nan_color
                elseif c < low && lowclip !== nothing
                    return lowclip
                elseif c > high && highclip !== nothing
                    return highclip
                else
                    Makie.interpolated_getindex(colormap, c, colorrange)
                end
            end
            return FaceIterator(cvec, faces)
        elseif color isa AbstractMatrix{<: Colorant} && uv !== nothing
            wsize = reverse(size(color))
            wh = wsize .- 1
            cvec = map(uv) do uv
                x, y = clamp.(round.(Int, Tuple(uv) .* wh) .+ 1, 1, wh)
                return color[end - (y - 1), x]
            end
            # TODO This is wrong and doesn't actually interpolate
            # Inside the triangle sampling the color image
            return FaceIterator(cvec, faces)
        end
    end
    error("Unsupported Color type: $(typeof(color))")
end

mesh_pattern_set_corner_color(pattern, id, c::Colorant) =
    Cairo.mesh_pattern_set_corner_color_rgba(pattern, id, rgbatuple(c)...)


text:
# Create file and initialises IOStream
fout = open("output.txt", "w")
a = write(fout, "First line of text.\n")
b = write(fout, "Second line of text.\n")
close(fout) # Closes stream, and writes to file.
println("$a, $b") # Write returns the number of characters written.

fout = open("integers.txt", "w")
num = 123
# The following line with not write the integer 123:
# write(fout, num)
# Convert to string before writing:
write(fout, string(num))
# Alternatively:
println(fout, string(456))
close(fout)

#=
Note the different I/O modes:
open(file, mode)
r: read
w: write, create, truncate
a: write, create. append
r+: read, write
w+: read, write, create, truncate
a+: read, write, create, append
=#

fout = open("output.txt", "a")
write(fout, "Third line of text.")
close(fout)

io = open("output.txt", "r")
println(read(io, String))
println(read(io, String))
println(read(io, String))
close(io)

io = IOBuffer();
write(io, "JuliaLang is a GitHub organization.")
truncate(io, 15)
println(String(take!(io)))
close(io)

# Paths and directories
println("homedir(): $(homedir())")
println("pwd(): $(pwd())")
println("readdir(): $(readdir())")
println("abspath(\".gitignore\"): $(abspath(".gitignore"))")
println("isfile(\".gitignore\"): $(isfile(".gitignore"))")
println("isdir(\".gitignore\"): $(isdir(".gitignore"))")
println("isdir(\"src/exercises\"): $(isdir("src/exercises"))")
println("ispath(\".gitignore\"): $(ispath(".gitignore"))")
println("ispath(\"src/exercises\"): $(ispath("src/exercises"))")

# Similar functions to unix utilities such as cd, rm, mkdir and mkpath. See https://docs.julialang.org/en/v1/base/file/#Base.Filesystem.walkdir
println("joinpath(\"/home/myuser\", \"example.jl\"): $(joinpath("/home/myuser", "example.jl"))")

function printallfiles(dir)
    for (root, dirs, files) in walkdir(dir)
        println("Directories in $root")
        for dir in dirs
            println(joinpath(root, dir)) # path to directories
        end
        println("Files in $root")
        for file in files
            println(joinpath(root, file)) # path to files
        end
    end
end

# printallfiles("src")

itr = walkdir("src");
println(first(itr))
println(first(itr))

# Remove files
rm("output.txt")
rm("integers.txt")


text:
function createabcfunction(prob, t, distancefunction, alg; kwargs...)
    function simfunc(params, constants, data)
        sol = solve(STANDARD_PROB_GENERATOR(prob, params), alg; saveat = t, kwargs...)
        if any((s.retcode != :Success for s in sol)) && any((s.retcode != :Terminated for s in sol))
            return Inf,nothing
        else
            simdata = convert(Array, sol)
            return distancefunction(data, simdata), nothing
        end
    end
end

function abc_inference(prob::DiffEqBase.DEProblem, alg, t, data, priors; ϵ=0.001,
                       distancefunction = euclidean, ABCalgorithm = ABCSMC, progress = false,
                       num_samples = 500, maxiterations = 10^5, kwargs...)

    abcsetup = ABCalgorithm(createabcfunction(prob, t, distancefunction, alg; kwargs...),
                            length(priors),
                            ϵ,
                            ApproxBayes.Prior(priors);
                            nparticles = num_samples,
                            maxiterations = maxiterations
                            )

    abcresult = runabc(abcsetup, data, progress = progress)
    return abcresult
end


text:
immutable BeliefDict <: Belief
  dict::Dict{String, Vector{Any}}
end

immutable BeliefSQLite <: Belief
  db::SQLite.DB
  table_name::String
end

get_beliefs(b::Belief, sample::String) = error("not implemented")

function get_beliefs(b::BeliefDict, sample::String; thresh=0)
  result = get(b.dict, sample, [])
  if thresh > 0
    result = filter(x->x[2]>thresh, result)
  end
  result
end

function get_beliefs(b::BeliefSQLite, sample::String; thresh=0)
  tbl = b.table_name
  query = "select term_name, prob from $tbl, terms, samples where $(tbl).term_id = terms.term_id " *
          " and $(tbl).sample_id = samples.sample_id and samples.sample_name = ? and prob>?"
  res = SQLite.query(b.db, query, values=[sample, thresh], nullable=false)
  collect(zip(res[:term_name], res[:prob]))
end

get_all_samples(b::Belief) = error("not implemented")
get_all_samples(b::BeliefDict) = keys(b.dict)

function get_all_samples(b::BeliefSQLite)
  tbl = b.table_name
  query = "select distinct sample_name from $tbl, samples where $(tbl).sample_id = samples.sample_id"
  res = SQLite.query(b.db, query, nullable=false)
  Set(res[:sample_name])
end

type BeliefSerializer
    db::SQLite.DB
    sample_dict::Dict{String, Int64}
    term_dict::Dict{String, Int64}
end

function BeliefSerializer(dbfn::String)
    db = SQLite.DB(dbfn)
    SQLite.execute!(db, "PRAGMA journal_mode=OFF;")
    SQLite.execute!(db, "PRAGMA mmap_size=268435456;")
    SQLite.execute!(db, "CREATE TABLE IF NOT EXISTS samples(sample_id INTEGER PRIMARY KEY, sample_name TEXT UNIQUE)")
    SQLite.execute!(db, "CREATE TABLE IF NOT EXISTS terms(term_id INTEGER PRIMARY KEY, term_name TEXT UNIQUE)")
    BeliefSerializer(db, load_dict(db, "samples"), load_dict(db, "terms"))
end

function load_dict(db, tblname)
    rows = SQLite.query(db, "select * from $tblname")
    Dict(zip(rows[:, 2], rows[:, 1]))
end

function get_id(val, db, dict, table_name, colname)
    val = string(val)
    get!(dict, val) do
      SQLite.query(db, "INSERT INTO $(table_name)($colname) VALUES (?)", values=[val])
      id = SQLite.query(db, "select last_insert_rowid()")[1][1]
      id
    end
end

get_sample_id(bs::BeliefSerializer, sample) = get_id(sample, bs.db, bs.sample_dict, "samples", "sample_name")
get_term_id(bs::BeliefSerializer, term) = get_id(term, bs.db, bs.term_dict, "terms", "term_name")

function store_beliefs(bs::BeliefSerializer, belief::BeliefSQLite, tblname; thresh=0.01, overwrite=false)
    if overwrite
      error("not supported")
    end
end
function store_beliefs(bs::BeliefSerializer, belief::BeliefDict, tblname; thresh=0.01, overwrite=true)
    has_table = tblname in SQLite.tables(bs.db)[:name]

    if has_table && !overwrite
      warn("requested no overwrite of $tblname which already exists")
      return
    end

    has_table && overwrite && SQLite.execute!(bs.db, "DROP TABLE $tblname")

    SQLite.execute!(bs.db, "CREATE TABLE $tblname (sample_id INTEGER, term_id INTEGER, prob REAL)")
    stmt = SQLite.Stmt(bs.db,  "INSERT INTO $tblname (sample_id, term_id, prob) VALUES (?, ?, ?)")
    for (sample, values) in belief.dict
        sid = get_sample_id(bs, sample)
        for (term, prob) in values
            prob < thresh && continue

            tid = get_term_id(bs, term)
            SQLite.bind!(stmt, [sid, tid, prob])
            SQLite.execute!(stmt)
        end
    end
    SQLite._close(stmt)
    SQLite.execute!(bs.db, "CREATE INDEX $(tblname)_sample_id_idx on $tblname (sample_id)")
    SQLite.execute!(bs.db, "CREATE VIEW $(tblname)_v(sample, term, prob) as select sample_name as sample, term_name as term, prob " *
                                        "from $(tblname), terms, samples where $(tblname).sample_id=samples.sample_id and " *
                                        "$(tblname).term_id=terms.term_id")
end

function store_params(db, params)
    params_str = JSON.json(params)
    SQLite.execute!(db, "CREATE TABLE IF NOT EXISTS meta(key text, value text)")
    SQLite.query(db, "INSERT INTO meta(key, value) VALUES (?, ?)", values=["params", params_str])
end

function store_beliefs(task::EZLearnTask, serializer::BeliefSerializer)
    gc() # otherwise transaction might fail due to stale statement objects
    SQLite.transaction(serializer.db) do
        store_params(serializer.db, task.params)

        for v in task.views
            id = v.id
            beliefs = task.beliefs[id]
            for (i, belief) in enumerate(beliefs)
                tbl_name = "$(id)_$i"
                println("Storing $(tbl_name)....")
                store_beliefs(serializer, belief, tbl_name)
                beliefs[i] = BeliefSQLite(serializer.db, tbl_name)
                gc()
                println("Done")
            end
        end
    end
end

store_beliefs(task::EZLearnTask, outfn::String) = store_beliefs(task, BeliefSerializer(outfn))

function intersect_labels_core(b1, b2, intersect_fn; keep_left=false, default_left=false, append_left=false, append_right=false)
    shared_samples = intersect(keys(b1), keys(b2))
    result = Dict{String, Vector{String}}()
    for sample in shared_samples
        terms = intersect_fn(b1[sample], b2[sample])
        if keep_left
          terms = union(terms, b1[sample])
        end
        if default_left && length(terms) == 0
          terms = b1[sample]
        end
        if length(terms) > 0
            result[sample] = terms
        end
    end
    if keep_left || default_left || append_left
      for sample in setdiff(keys(b1), keys(b2))
        result[sample] = b1[sample]
      end
    end
    if append_right
      for sample in setdiff(keys(b2), keys(b1))
        result[sample] = b2[sample]
      end
    end
    result
end

intersect_labels_simple(b1, b2) = intersect_labels_core(b1, b2, intersect)

OBOParse.descendants(o::Ontology, terms::Vector{Term}, syms::Vector{Symbol}) = union([descendants(o, t, syms) for t in terms]...)

function ontology_intersector(ontology, symbols)
    function f(vals1, vals2)
        terms1 = Term[gettermbyid(ontology, t) for t in vals1]
        terms2 = Term[gettermbyid(ontology, t) for t in vals2]

        terms1_all = union(terms1, descendants(ontology, terms1, symbols))
        terms2_all = union(terms2, descendants(ontology, terms2, symbols))

        resvals = union(intersect(terms1, terms2_all), intersect(terms2, terms1_all))

        return String[t.id for t in resvals]
    end
    return f
end

function threshold_beliefs{T<:Belief}(beliefs::T, threshold)
    result = Dict{String, Vector{String}}()
    for sample in get_all_samples(beliefs)
        val = get_beliefs(beliefs, sample; thresh=threshold)
        if length(val) > 0
          result[sample] = String[t[1] for t in val]
        end
    end
    result
end

function remove_redundant(ontology, terms::Vector{Term}, rels::Vector{Symbol})
    if length(terms) <= 1
        return terms
    end

    terms = unique(terms)
    termSet = Set(terms)
    n = length(terms)
    keep = ones(Bool, n)
    for i in 1:n
        desc = descendants(ontology, terms[i], rels)
        keep[i] = length(intersect(desc, termSet)) == 0
    end

    terms[keep]
end

getterm(t::Integer, ontology) = gettermbyid(ontology, t)
function getterm(t::String, ontology)
    if haskey(ontology.terms, t)
        return ontology[t]
    else
        return gettermbyid(ontology, parse(Int, t))
    end
end

function threshold_beliefs_nonred(ontology, beliefs, threshold, rels)
    result = Dict{String, Vector{String}}()
    for gsm in get_all_samples(beliefs)
        terms = get_beliefs(beliefs, gsm; thresh=threshold)
        if length(terms) > 0
            #termsids = String[t[1] for t in terms]
            termobjs = Term[getterm(t[1], ontology) for t in terms]
            nonred = remove_redundant(ontology, termobjs, rels)
            @assert length(nonred) > 0
            result[gsm] = String[t.id for t in nonred]
        end
    end
    result
end

function remove_redundant_beliefs(ontology, beliefs, rels)
    result = Dict{String, Vector{String}}()
    for gsm in keys(beliefs)
        terms = beliefs[gsm]
        if length(terms) > 0
            termobjs = Term[gettermbyid(ontology, t) for t in terms]
            nonred = remove_redundant(ontology, termobjs, rels)
            result[gsm] = String[t.id for t in nonred]
        end
    end
    result
end


text:
struct LaplaceApproximation{SM, F,O,YT,XT,WT}
    size::Tuple{Int, Int}
    Q::SM # Sparse precision matrix
    # H and chol in obj
    # H::SM # Hessian
    # chol::CH # Symbolic factorization of Q
    lik::F # loglikelihood, derivative and second derivative.
    obj::O # Optim.TwiceDifferentiable
    y::YT # Response
    x::XT # Fixed effect ???
    w::WT # Integration weights
end
function poissonlikelihood(y, x, z)
    exz = exp(x+z)
    y*(x+z) - exz, y - exz, -exz
end

function poisson_gp(locationweight, fixed_field, prec, range, nu, integration_weights, gridsize)
    Q = Qmatrix_matern2d(gridsize, prec, range, nu)
    chol = cholesky(Q)
    Qlogdet = logdet(chol)

    la = LaplaceApproximation(gridsize, locationweight, fill(integration_weights, gridsize),
        Q.data, chol)
    poisson_gp(la, fixed_field, Q, Qlogdet)
end
function LA_factors(zhat, H, Q, Qlogdet)
    d = length(zhat)

    pz = 0.5*Qlogdet - d/2*log(2*pi) - 0.5*(zhat'*Q*zhat)

    m = d/2*log(2*pi) - 0.5*logdet(H)

    a = m + pz
    ifelse(isnan(a), -Inf, a), m, pz
end

function poisson_gp(la::LaplaceApproximation, fixed_field, Q, Qlogdet)
    zhat, H = poisson_hatZ(la, fixed_field, Q)
    # H is the hessian of the negative loglikelihood
    x = zhat

    l = LA_factors(zhat, H, Q, Qlogdet)
    l, reshape(zhat, la.size), H, Q
end
function LaplaceApproximation(gridsize, locationweight, integration_weights, Q, chol)
    H = copy(Q)
    y = locationweight
    x = zeros(gridsize)
    w = integration_weights

    obj = LA_poisson_objective_chol(y, x, w, Q, H, chol)
    LaplaceApproximation(gridsize, Q, #copy(Q), cholesky(Q),
        nothing, obj, y, x, w)
end

function poisson_hatZ(la::LaplaceApproximation, fixed_field::Number, Q)
    poisson_hatZ(la, FillArrays.Fill(fixed_field, la.size), Q)
end
function poisson_hatZ(la::LaplaceApproximation, fixed_field, Q)
    gridsize = size(fixed_field)
    copy!(la.x, fixed_field)

    copy!(la.Q, Q.data)
    initial = vec(initial_guess(la.x, la.y, Q[1]))
    obj = la.obj

    o = try
        Optim.optimize(obj, initial)
    catch e
        if e isa AssertionError
            println(e)
        end
        rethrow()
    end
    if !Optim.converged(o)
        @warn("Newton Not converged.")
    end
    x = Optim.minimizer(o)
    H = Optim.hessian!(obj, x)

    x, H
end

# Here Q is sparse matrix
function LA_poisson_objective_chol(locationweight, fixed_field, integration_weights, Q, H, cholQ)
    function fgh!(F, D, chol, x)
        # _negative_ log likelihood
        random_field = x
        if chol !== nothing
            copy!(H, Q)
        end
        if D !== nothing || F !== nothing
            Qx = Q*x
            if D !== nothing
                D .= Qx
                weightstonodes!(D, locationweight, -1)
            end
            if F !== nothing
                xy = sum(locationweight) do (i, w)
                    w*random_field[i]
                end
                F = 0.5 * x'*Qx - xy
            end
        end
        for i in eachindex(random_field)
            a = integration_weights[i] * exp(fixed_field[i] + random_field[i])
            if !isfinite(a)

                @warn("Infinite a")
                return -Inf
            end
            F !== nothing && (F += a)
            D !== nothing && (D[i] += a)
            chol !== nothing && (H[i, i] += a)
        end
        if chol !== nothing
            cholesky!(chol.data, Symmetric(H))
        end
        if F !== nothing
            return F
        end
    end
    initial = vec(zero(fixed_field))
    Optim.TwiceDifferentiable(Optim.only_fgh!(fgh!), initial, 0.0, copy(initial), OptimCholWrapper(cholQ))
end

# If Q is a dense matrix some things are different
function LA_poisson_objective_chol(locationweight, fixed_field, integration_weights, Q::Matrix{Float64}, H, cholQ)
    function fgh!(F, D, H, x)
        # _negative_ log likelihood
        random_field = x
        if H !== nothing
            copy!(H, Q)
        end
        if D !== nothing || F !== nothing
            Qx = Q*x
            if D !== nothing
                D .= Qx
                weightstonodes!(D, locationweight, -1)
            end
            if F !== nothing
                xy = sum(locationweight) do (i, w)
                    w*random_field[i]
                end
                F = 0.5 * x'*Qx - xy
            end
        end
        for i in eachindex(random_field)
            a = integration_weights[i] * exp(fixed_field[i] + random_field[i])
            if !isfinite(a)

                @warn("Infinite a")
                return -Inf
            end
            F !== nothing && (F += a)
            D !== nothing && (D[i] += a)
            H !== nothing && (H[i, i] += a)
        end
        return F
    end
    initial = vec(zero(fixed_field))
    Optim.TwiceDifferentiable(Optim.only_fgh!(fgh!), initial, 0.0, copy(initial), copy(Q))
end


text:
module DemoPackageJCM

print_greeting(io::IO=stdout)=printstyled("Hello, moon!"; color =:green)

end


text:
function creategrid(;xrange::(Real, Real)=(-1,1), yrange::(Real, Real)=(-1,1), zrange::(Real, Real)=(-1,1), camera=pcamera)
	xyplane = genquad(Vec3(0, 0, 0), Vec3(1, 0, 0), Vec3(0, 1, 0))
	zyplane = genquad(Vec3(0, 0, 0), Vec3(0, 0, 1.2), Vec3(0, 1, 0))
	zxplane = genquad(Vec3(0, 1, 0), Vec3(0, 0, 1.2), Vec3(1, 0, 0))

	v,uv,n,i = mergemesh(xyplane, zyplane, zxplane)

	grid = RenderObject(@compat(Dict(
			:vertexes 			  	=> GLBuffer(v),
			:indexes			   	=> indexbuffer(i),
			#:grid_color 		  => Float32[0.1,.1,.1, 1.0],
			:bg_color 			  	=> Input(Vec4(1, 1, 1, 0.01)),
			:grid_thickness  		=> Input(Vec3(2)),
			:gridsteps  		  	=> Input(Vec3(10)),
			:mvp 				    => camera.projectionview
		)), gridshader)
	prerender!(grid, glDisable, GL_DEPTH_TEST, glDisable, GL_CULL_FACE, enabletransparency)
	postrender!(grid, render, grid.vertexarray)
	return grid
end

initgrid() = global gridshader = TemplateProgram(joinpath(shaderdir,"grid.vert"), joinpath(shaderdir,"grid.frag"))
init_after_context_creation(initgrid)

export creategrid


text:
# This file was generated, do not modify it. # hide
using MLJ
using  PrettyPrinting
import DataFrames
import Statistics
MLJ.color_off() # hide

X, y = @load_reduced_ames
X = DataFrames.DataFrame(X)
@show size(X)
first(X, 3) |> pretty


text:
using Latexify
using Test
import Latexify: PlainNumberFormatter, PrintfNumberFormatter

@test FancyNumberFormatter() == FancyNumberFormatter(4) == FancyNumberFormatter("%.4g", "\\cdot") == FancyNumberFormatter("%.4g", s"\g<mantissa> \\cdot 10^{\g<exp>}")

x = -23.4728979e7

@test PlainNumberFormatter()(x) == "-2.34728979e8"
@test PrintfNumberFormatter("%.4g")(x) == "-2.347e+08"
@test StyledNumberFormatter()(x) == "-2.347 \\mathrm{e} 8"
@test FancyNumberFormatter()(x) == "-2.347 \\cdot 10^{8}"
@test FancyNumberFormatter("%.5E", s"\g<mantissa>,\g<before_dp>,\g<sign>,\g<before_dp_nosign>,\g<after_dp>,\g<e_or_E>,\g<raw_exp>,\g<exp>,\g<sign_exp>")(x) == "-2.34729,-2,-,2,34729,E,+08,8,"

y = 0xf43

@test StyledNumberFormatter()(y) == FancyNumberFormatter()(y) == "\\mathtt{0x0f43}"


text:
for F in (:diag, :diagm, :logdet, :norm, :normalize, :normalize!, :pinv, :qr, :transpose)
   @eval $F(x::A) where {T<:AbstractFinite, A<:AbstractVector{T}} = T.($F(float(T).(x)))
end

for F in (:adjoint, :diag, :diagm, :expm, :exmp!, :logdet, :norm, :normalize, :normalize!,
          :nullspace, :pinv, :qr, :transpose)
   @eval $F(x::V) where {T<:AbstractFinite, V<:Vector{T}} = T.($F(float(T).(x)))
end

for F in (:adjoint, :bunchkaufman, :cond, :condskeel, :det, :diag, :diagind,
          :eigmax, :eigmin, :eigvecs, :expm, :expm!,
          :isdiag, :ishermitian, :isposdef, :isposdef!, :issymmetric, :istril, :istriu,
          :logabsdet, :logdet, :lu, :norm, :opnorm, :qr, :rank, :svdvals,
          :tr, :transpose, :tril, :tril!, :triu, :triu!)
   @eval $F(x::M) where {T<:AbstractFinite, M<:AbstractMatrix{T}} = T.($F(float(T).(x)))
end

for F in (:bunchkaufman, :cholesky, :cholesky!, :cond, :condskeel,
          :det, :diag, :diagind, :eigen, :eigmax, :eigmin, :eigvals, :eigvecs,
          :expm, :expm!, :factorize, :hessenberg, :isdiag, :ishermitian,
          :isposdef, :isposdef!, :issymmetric, :istril, :istriu, :logabsdet,
          :logdet, :lu, :lu!, :norm, :nullspace, :opnorm, :pinv, :qr, :qr!,
          :rank, :schur, :svd, :svdvals, :tr, :transpose, :tril, :tril!, :triu, :triu!)
   @eval $F(x::M) where {T<:AbstractFinite, M<:Matrix{T}} = T.($F(float(T).(x)))
end

for F in (:adjoint, :expm, :exmp!, :logdet, :norm, :transpose)
   @eval $F(x::A) where {T<:AbstractFinite, A<:AbstractArray{T}} = T.($F(float(T).(x)))
end

for F in (:adjoint!, :copyto!, :dot, :eigen, :eigvals, :eigvecs, :expm, :expm!,
          :fill!, :kron, :qr, :schur, :svd, :svdvals, :transpose!)
    @eval $F(x::M, y::M) where {T<:AbstractFinite, M<:AbstractMatrix{T}} = T.($F(float(T).(x), float(T).(y)))
end

#=
stack overflow if used

for F in (:copyto!, :cross, :dot, :expm, :expm!, :fill!, :kron)
    @eval $F(x::V, y::V) where {T<:AbstractFinite, V<:AbstractVector{T}} = T.($F(float(T).(x), float(T).(y)))
end

=#


text:
using SafeTestsets

@safetestset "add_nums test sets" begin
  include("add_nums.jl")
end

@safetestset "cat_str test sets" begin
  include("cat_str.jl")
end


text:
using ModelingToolkit, IfElse, LinearAlgebra, StaticArrays
import LinearAlgebra.norm
norm(x::AbstractVector{Num}) = sqrt(sum(x.^2))

# attempt to create symbolic version of GeometryPrimitives.Polygon

# nv vertices in 2D plane
nv = 4  # first try trapazoid like in ridge waveguide
@parameters w, t, θ # compare with wₜₒₚ,t_core,θ
# derived parameters for vertices
wt_half = w / 2
wb_half = wt_half + ( t * tan(θ) )
tc_half = t / 2
# vertices
v =    	[   wt_half     tc_half
            -wt_half    tc_half
            -wb_half    -tc_half
            wb_half     -tc_half    ]
vs = SMatrix{4,2}(v)
Δv = v - circshift(v,1)
Δvs = SMatrix{4,2}(Δv)

one_mone = [1, -1]
one_mones = SVector(1, -1)
@tullio n[i,j] := Δv[i,3-j] * $one_mone[j] / sqrt(Δv[i,1]^2 + Δv[i,2]^2 ) nograd=one_mone
ns = SMatrix{4,2}(n)
lbnd = transpose(minimum(v;dims=1))
ubnd = transpose(maximum(v;dims=1))

lbnds = transpose(minimum(vs;dims=1))
ubnds = transpose(maximum(vs;dims=1))
bnds = lbnd, ubnd

function Δxe(x,v,n)
    @tullio res[k] := n[k,a] * ( x[a] - v[k,a] )
end

minbool(x) = map(a->(a <= minimum(x)), x)
maxbool(x) = map(a->(a >= maximum(x)), x)
minone(x) = IfElse.ifelse.(x .<= minimum(x), 1, 0)
maxone(x) = IfElse.ifelse.(x .>= maximum(x), 1, 0)
atmin(x,y) = sum(IfElse.ifelse.(y .<= minimum(y), x, 0))
atmax(x,y) = sum(IfElse.ifelse.(y .>= maximum(y), x, 0))

# test atmin, atmax
# @parameters u[1:10], q[1:10]
# fminbool,fminbool! = eval.(build_function(minbool(q),q))
# fmaxbool,fmaxbool! = eval.(build_function(maxbool(q),q))
# fatmin = eval(build_function(atmin(u,q),u,q))
# fatmax = eval(build_function(atmax(u,q),u,q))
# un = [1,2,3,4,5,6,7,8,9,10]
# qn = [3,4,3,2,10,5,3,5,1,2]
# fmaxbool(un)
# fatmin(un,qn)
# fatmax(un,qn)

function pgon_cout2(x,v,n)
	Δxv = x' .- v
	lΔxv = sqrt.(sum(x->x^2,dxv,dims=2))
	imin = minone(lΔxv)
	surf = sum(v.*imin,dims=1)[1,:]
	# IfElse.ifelse(obd(x)&&)
	# imin₋₁ = circshift(imin,-1) #for ifelse if needed
	n0 = x-surf
	# return surf, n0./norm(n0)
	hcat(surf, n0./norm(n0))
end

function pgon_cout_not2(x,v,n,Δxe)
	imax = maxone(Δxe)
	vmax = sum(imax.*v,dims=1)[1,:]
	nout = sum(imax.*n,dims=1)[1,:]
	Δx = (nout⋅(vmax-x)) .* nout
    surf = x + Δx
	# return surf, nout
	hcat(surf, nout)
end

function spn_pgon(x,v,n)
    # @tullio Δxe[k] := n[k,a] * ( x[a] - v[k,a] )
	Δxe = sum(n .* (x' .- v), dims=2)[:,1]
    # cout = sum( 1 .- signbit.( Δxe ) )
	# cout = sum( Δxe .> 0 )
	cout = sum(IfElse.ifelse.(Δxe .> 0, 1, 0))
	res = IfElse.ifelse(
		isequal(cout,2),
		pgon_cout2(x,v,n),
		pgon_cout_not2(x,v,n,Δxe)
		)
	return res #::Matrix
end

# fspn_pgon,fspn_pgon! = eval.(build_function(spn_pgon(x,v,n),x,w,t,θ))
fspn_pgon,fspn_pgon! = eval.(build_function(spn_pgon(x,v,n),x,ps))
fDspn_pgon = eval(build_function(Dspn_pgon(x,v,n),x,ps))
fspn_pgons,fspn_pgons! = eval.(build_function(spn_pgon(xs,vs,ns),x,w,t,θ))

fpgon_cout2,fpgon_cout2! = eval.(build_function(pgon_cout2(x,v,n),x,w,t,θ))
fpgon_cout_not2,fpgon_cout_not2! = eval.(build_function(pgon_cout_not2(x,v,n,Δxe(x,v,n)),x,w,t,θ))

xnum = [0.3,0.5]
@parameters x[1:2]
spn_pgon(xnum,v,n)
spn_pgon(x,v,n)

f1 = eval(build_function(spn_pgon(xnum,v,n),w,t,θ))
f2,f2! = eval.(build_function(Δxe(xnum,v,n),w,t,θ))
f3,f3! = eval.(build_function(Δxe(x,v,n),x,w,t,θ))
f4 = eval(build_function(spn_pgon(x,v,n),x,w,t,θ))
flbnd,flbnd! = eval.(build_function(lbnd,w,t,θ))
fubnd,fubnd! = eval.(build_function(ubnd,w,t,θ))
fbnd(a,b,c) = flbnd(a,b,c),fubnd(a,b,c)
# @variables v[nv,2]

# compare with original functions
using GeometryPrimitives, StaticArrays, LinearAlgebra
function pg(w,t,θ)
    wt_half = w / 2
    wb_half = wt_half + ( t * tan(θ) )
    tc_half = t / 2
    # vertices
    verts = SMatrix{4,2}( [  wt_half     tc_half
                            -wt_half    tc_half
                            -wb_half    -tc_half
                            wb_half     -tc_half    ] )
    GeometryPrimitives.Polygon(					                        # Instantiate 2D polygon, here a trapazoid
                    verts,			                            # v: polygon vertices in counter-clockwise order
                    2.4,					                                    # data: any type, data associated with box shape
                )
end

bounds(pg1)
function _Δxe_poly(x::SVector{2},v::SMatrix{K,2},n::SMatrix{K,2})::SVector{K}  where {K} #,T<:Real}
	@tullio out[k] := n[k,a] * ( x[a] - v[k,a] ) # edge line eqns for a K-point Polygon{K} `s`
end
function Δxe_poly(x::AbstractVector{<:Real},s::Polygon{K})::SVector{K} where K
	_Δxe_poly(x,s.v,s.n)
end
function onbnd(Δxe::SVector{K},s::Polygon{K})::SVector{K,Bool} where K
	map(x->abs(x)≤s.rbnd,Δxe)
end
function onbnd(x::SVector{2},s::Polygon)::SVector{K,Bool} where K
	map(x->abs(x)≤s.rbnd,Δxe_poly(x,s))
end
function cout(x::SVector{2},s::Polygon)::Int
	mapreduce((a,b)->Int(a|b),+,onbnd(x,s),map(isposdef,Δxe_poly(x,s)))
end
function cout(Δxe::SVector{K,<:Real},obd::SVector{K,Bool})::Int where K
	mapreduce((a,b)->Int(a|b),+,obd,map(isposdef,Δxe))
end

pg1 = pg(1.7,0.7,π/14.0)
# x = SVector(0.2,0.3)
∆xe1 = _Δxe_poly(SVector{2}(xnum),pg1.v,pg1.n)
obd1 = onbnd(∆xe1,pg1)
co1 = cout(∆xe1,obd1)

## reference code

p = [
       1.7,                #   top ridge width         `w_top`         [μm]
       0.7,                #   ridge thickness         `t_core`        [μm]
       π / 14.0,           #   ridge sidewall angle    `θ`             [radian]
       2.4,                #   core index              `n_core`        [1]
       1.4,                #   substrate index         `n_subs`        [1]
       0.5,                #   vacuum gap at boundaries `edge_gap`     [μm]
               ];
Δx,Δy,Δz,Nx,Ny,Nz = 6.0, 4.0, 1.0, 128, 128, 1;
rwg(p) = ridge_wg(p[1],p[2],p[3],p[6],p[4],p[5],Δx,Δy)
ms = ModeSolver(1.45, rwg(p), Δx, Δy, Δz, Nx, Ny, Nz)

function ridge_wg2(wₜₒₚ::T1,t_core::T1,θ::T1,edge_gap::T1,ε_core,ε_subs,Δx::T2,Δy::T2)::Vector{<:GeometryPrimitives.Shape} where {T1<:Real,T2<:Real}
    t_subs = (Δy -t_core - edge_gap )/2.
    c_subs_y = -Δy/2. + edge_gap/2. + t_subs/2.
    # ε_core = ε_tensor(n_core)
    # ε_subs = ε_tensor(n_subs)
    wt_half = wₜₒₚ / 2
    wb_half = wt_half + ( t_core * tan(θ) )
    tc_half = t_core / 2
    verts =     [   wt_half     -wt_half     -wb_half    wb_half
                    tc_half     tc_half    -tc_half      -tc_half    ]'
    core = GeometryPrimitives.Polygon(					                        # Instantiate 2D polygon, here a trapazoid
                    SMatrix{4,2}(verts),			                            # v: polygon vertices in counter-clockwise order
                    ε_core,					                                    # data: any type, data associated with box shape
                )
    ax = [      1.     0.
                0.     1.      ]
    b_subs = GeometryPrimitives.Box(					                # Instantiate N-D box, here N=2 (rectangle)
                    [0. , c_subs_y],            	# c: center
                    [Δx - edge_gap, t_subs ],	# r: "radii" (half span of each axis)
                    ax,	    		        # axes: box axes
                    ε_subs,					        # data: any type, data associated with box shape
                )
    return [core,b_subs]
end

rwg2(p) = ridge_wg2(p[1],p[2],p[3],p[6],MgO_LiNbO₃,SiO₂,Δx,Δy)
shapes = rwg(p)
shapes2 = rwg2(p)


text:
using AccuracyAtTop
using EvalMetrics
using Flux
using MLDatasets
using Plots
using ProgressMeter

using Flux: gpu

# ------------------------------------------------------------------------------------------
# Auxiliary functions
# ------------------------------------------------------------------------------------------
function reshape_data(X::AbstractArray{T, 3}, y::AbstractVector, pos_class) where T
    s = size(X)
    return reshape(X, s[1], s[2], 1, s[3]), reshape(y, 1, :) .== pos_class
end

function eval_scores(model, batches)
    targets = Bool[]
    scores = Float32[]

    @showprogress for (x,y) in batches
        append!(targets, vec(y))
        append!(scores, vec(model(x)))
    end
    return targets, scores
end

# ------------------------------------------------------------------------------------------
# Data preparation
# ------------------------------------------------------------------------------------------
using Flux.Data: DataLoader

T = Float32
pos_class = 0
batchsize = 128
device = gpu

X_train, y_train = reshape_data(MLDatasets.FashionMNIST.traindata(T)..., pos_class);
X_test, y_test = reshape_data(MLDatasets.FashionMNIST.testdata(T)..., pos_class);

batches_train = (device(batch) for batch in DataLoader((X_train, y_train); batchsize))
batches_test = (device(batch) for batch in DataLoader((X_test, y_test); batchsize))

# ------------------------------------------------------------------------------------------
# Model preparation
# ------------------------------------------------------------------------------------------
model = Chain(
    # First convolution
    Conv((5, 5), 1=>20, stride=(1,1), relu),
    MaxPool((2,2)),

    # Second convolution
    Conv((5, 5), 20=>50, stride=(1,1), relu),
    MaxPool((2,2)),

    flatten,
    Dense(4*4*50, 500),
    Dense(500, 1)
) |> device

# objective
surrogate = quadratic;
pars = params(model);
aatp = DeepTopPush()

loss(x, y) = objective(aatp, y, model(x); surrogate)

# ------------------------------------------------------------------------------------------
# Model training
# ------------------------------------------------------------------------------------------
opt = ADAM(0.001);
n_epochs = 10

@showprogress for i in 1:n_epochs
    Flux.train!(loss, pars, batches_train, opt)
end

# ------------------------------------------------------------------------------------------
# Evaluation
# ------------------------------------------------------------------------------------------
tar_train, s_train = eval_scores(model, batches_train)
tar_test, s_test = eval_scores(model, batches_test)

kwargs = (xscale = :log10, xlims = (1e-5, 1))
plot(
    rocplot(tar_train, s_train; title = "Train", kwargs...),
    rocplot(tar_test, s_test; title = "Test", kwargs...),
    size = (800, 400)
)


text:
using LiveServer
servedocs(doc_env=true, foldername=@__DIR__)


text:
function f!(resid, x)
    for i in eachindex(x)
       resid[i] = sin(x[i]) + x[i]^3
    end
end
x = ones(5)
@test Sundials.kinsol(f!, x, linear_solver=:Band, jac_upper=0, jac_lower=0) == Sundials.kinsol(f!, x)


text:
function number_of_trees(right,down)

   data= readdlm("input_day3", '\n')
   horizontal = div( (size(data)[1])*right, down )
   lline = length(data[1])

   #number of times to repet the pattern
   rep = div(horizontal,lline)+1

   #repeat the pattern the ammount needed
   for i in eachindex(data)
      data[i] = data[i]^rep
   end

   #horizontal index
   h_i = 1

   #number of trees
   n_t = 0

   for i in Iterators.countfrom(1,down)
      i > size(data)[1]-down && break
      line = data[i+down]
      h_i += right
      index = h_i

      if line[h_i] == '#'
        n_t += 1
      end
   end

   return (n_t)

end

using DelimitedFiles

r = [ 1, 3, 5, 7, 1 ]
d = [ 1, 1, 1, 1, 2 ]
slope = zeros(5)

product = 1

for i in eachindex(r)
   slope[i] = number_of_trees(r[i], d[i])
   global product *= slope[i]
   println(product)
end

println(product)