text:
struct Fred
    a::Float64
end
Base.zero(::Fred) = Fred(0.0)
Base.zero(::Type{Fred}) = Fred(0.0)

struct Freddy{T,N}
    a::Array{T,N}
end
Base.:(==)(a::Freddy, b::Freddy) = a.a == b.a

struct Mary
    a::Fred
end

struct TwoFields
    a::Float64
    b::Float64
end

@testset "projection" begin
    @testset "display" begin
        @test startswith(repr(ProjectTo(Fred(1.1))), "ProjectTo{Fred}(")
        @test repr(ProjectTo(1.1)) == "ProjectTo{Float64}()"
    end

    @testset "fallback" begin
        @test Fred(1.2) == ProjectTo(Fred(1.1))(Fred(1.2))
        @test Fred(0.0) == ProjectTo(Fred(1.1))(ZeroTangent())
        @test Fred(3.2) == ProjectTo(Fred(1.1))(@thunk(Fred(3.2)))
        @test Fred(1.2) == ProjectTo(Fred(1.1))(Tangent{Fred}(; a=1.2))

        # struct with complicated field
        x = Freddy(zeros(2, 2))
        dx = Tangent{Freddy}(; a=ZeroTangent())
        @test x == ProjectTo(x)(dx)

        # nested structs
        f = Fred(0.0)
        tf = Tangent{Fred}(; a=ZeroTangent())
        m = Mary(f)
        dm = Tangent{Mary}(; a=tf)
        @test m == ProjectTo(m)(dm)

        # two fields
        tfa = TwoFields(3.0, 0.0)
        tfb = TwoFields(0.0, 3.0)
        @test tfa == ProjectTo(tfa)(Tangent{TwoFields}(; a=3.0))
        @test tfb == ProjectTo(tfb)(Tangent{TwoFields}(; b=3.0))
    end

    @testset "to Real" begin
        # Float64
        @test 3.2 == ProjectTo(1.0)(3.2)
        @test 0.0 == ProjectTo(1.0)(ZeroTangent())
        @test 3.2 == ProjectTo(1.0)(@thunk(3.2))

        # down
        @test 3.2 == ProjectTo(1.0)(3.2 + 3im)
        @test 3.2f0 == ProjectTo(1.0f0)(3.2)
        @test 3.2f0 == ProjectTo(1.0f0)(3.2 - 3im)

        # up
        @test 2.0 == ProjectTo(1.0)(2.0f0)
    end

    @testset "to Number" begin
        # Complex
        @test 2.0 + 4.0im == ProjectTo(1.0im)(2.0 + 4.0im)

        # down
        @test 2.0 + 0.0im == ProjectTo(1.0im)(2.0)
        @test 0.0 + 0.0im == ProjectTo(1.0im)(ZeroTangent())
        @test 0.0 + 0.0im == ProjectTo(1.0im)(@thunk(ZeroTangent()))

        # up
        @test 2.0 + 0.0im == ProjectTo(1.0im)(2.0)
    end

    @testset "to Array" begin
        # to an array of numbers
        x = zeros(2, 2)
        @test [1.0 2.0; 3.0 4.0] == ProjectTo(x)([1.0 2.0; 3.0 4.0])
        @test x == ProjectTo(x)(ZeroTangent())

        x = zeros(2)
        @test x == ProjectTo(x)(@thunk(ZeroTangent()))

        x = zeros(Float32, 2, 2)
        @test x == ProjectTo(x)([0.0 0; 0 0])

        x = [1.0 0; 0 4]
        @test x == ProjectTo(x)(Diagonal([1.0, 4]))

        # to a array of structs
        x = [Fred(0.0), Fred(0.0)]
        @test x == ProjectTo(x)([Fred(0.0), Fred(0.0)])
        @test x == ProjectTo(x)([ZeroTangent(), ZeroTangent()])
        @test x == ProjectTo(x)([ZeroTangent(), @thunk(Fred(0.0))])
        @test x == ProjectTo(x)(ZeroTangent())
        @test x == ProjectTo(x)(@thunk(ZeroTangent()))

        x = [Fred(1.0) Fred(0.0); Fred(0.0) Fred(4.0)]
        @test x == ProjectTo(x)(Diagonal([Fred(1.0), Fred(4.0)]))
    end

    @testset "To Array of Arrays" begin
        # inner arrays have same type but different sizes
        x = [[1.0, 2.0, 3.0], [4.0, 5.0]]
        @test x == ProjectTo(x)(x)
        @test x == ProjectTo(x)([[1.0 + 2im, 2.0, 3.0], [4.0 + 2im, 5.0]])

        # This makes sure we don't fall for https://github.com/JuliaLang/julia/issues/38064
        @test [[0.0, 0.0, 0.0], [0.0, 0.0]] == ProjectTo(x)(ZeroTangent())
    end

    @testset "Array{Any} with really messy contents" begin
        # inner arrays have same type but different sizes
        x = [[1.0, 2.0, 3.0], [4.0 + im 5.0], [[[Fred(1)]]]]
        @test x == ProjectTo(x)(x)
        @test x == ProjectTo(x)([[1.0 + im, 2.0, 3.0], [4.0 + im 5.0], [[[Fred(1)]]]])
        # using a different type for the 2nd element (Adjoint)
        @test x == ProjectTo(x)([[1.0 + im, 2.0, 3.0], [4.0 - im, 5.0]', [[[Fred(1)]]]])

        @test [[0.0, 0.0, 0.0], [0.0im 0.0], [[[Fred(0)]]]] == ProjectTo(x)(ZeroTangent())
    end

    @testset "to Diagonal" begin
        d_F64 = Diagonal([0.0, 0.0])
        d_F32 = Diagonal([0.0f0, 0.0f0])
        d_C64 = Diagonal([0.0 + 0im, 0.0])
        d_Fred = Diagonal([Fred(0.0), Fred(0.0)])

        # from Matrix
        @test d_F64 == ProjectTo(d_F64)(zeros(2, 2))
        @test d_F64 == ProjectTo(d_F64)(zeros(Float32, 2, 2))
        @test d_F64 == ProjectTo(d_F64)(zeros(ComplexF64, 2, 2))

        # from Diagonal of Numbers
        @test d_F64 == ProjectTo(d_F64)(d_F64)
        @test d_F64 == ProjectTo(d_F64)(d_F32)
        @test d_F64 == ProjectTo(d_F64)(d_C64)

        # from Diagonal of AbstractTangent
        @test d_F64 == ProjectTo(d_F64)(ZeroTangent())
        @test d_C64 == ProjectTo(d_C64)(ZeroTangent())
        @test d_F64 == ProjectTo(d_F64)(@thunk(ZeroTangent()))
        @test d_F64 == ProjectTo(d_F64)(Diagonal([ZeroTangent(), ZeroTangent()]))
        @test d_F64 == ProjectTo(d_F64)(Diagonal([ZeroTangent(), @thunk(ZeroTangent())]))

        # from Diagonal of structs
        @test d_Fred == ProjectTo(d_Fred)(ZeroTangent())
        @test d_Fred == ProjectTo(d_Fred)(@thunk(ZeroTangent()))
        @test d_Fred == ProjectTo(d_Fred)(Diagonal([ZeroTangent(), ZeroTangent()]))

        # from Tangent
        @test d_F64 == ProjectTo(d_F64)(Tangent{Diagonal}(; diag=[0.0, 0.0]))
        @test d_F64 == ProjectTo(d_F64)(Tangent{Diagonal}(; diag=[0.0f0, 0.0f0]))
        @test d_F64 == ProjectTo(d_F64)(
            Tangent{Diagonal}(; diag=[ZeroTangent(), @thunk(ZeroTangent())])
        )
    end

    @testset "to $SymHerm" for SymHerm in (Symmetric, Hermitian)
        data = [1.0+1im 2-2im; 3 4]

        x = SymHerm(data)
        @test x == ProjectTo(x)(data)
        @test x == ProjectTo(x)(Tangent{typeof(x)}(; data=data, uplo=NoTangent()))

        x = SymHerm(data, :L)
        @test x == ProjectTo(x)(data)

        data = [1.0-2im 0; 0 4]
        x = SymHerm(data)
        @test x == ProjectTo(x)(Diagonal([1.0 - 2im, 4.0]))

        data = [0.0+0im 0; 0 0]
        x = SymHerm(data)
        @test x == ProjectTo(x)(ZeroTangent())
        @test x == ProjectTo(x)(@thunk(ZeroTangent()))
    end

    @testset "to $UL" for UL in (UpperTriangular, LowerTriangular)
        data = [1.0+1im 2-2im; 3 4]

        x = UL(data)
        @test x == ProjectTo(x)(data)
        @test x == ProjectTo(x)(Tangent{typeof(x)}(; data=data))

        data = [0.0+0im 0; 0 0]
        x = UL(data)
        @test x == ProjectTo(x)(Diagonal(zeros(2)))
        @test x == ProjectTo(x)(ZeroTangent())
        @test x == ProjectTo(x)(@thunk(ZeroTangent()))
    end

    @testset "to Transpose" begin
        x = rand(ComplexF64, 3, 4)
        t = transpose(x)
        mt = collect(t)
        a = adjoint(x)
        ma = collect(a)

        @test t == ProjectTo(t)(mt)
        @test conj(t) == ProjectTo(t)(ma)
        @test zeros(4, 3) == ProjectTo(t)(ZeroTangent())
        @test zeros(4, 3) == ProjectTo(t)(Tangent{Transpose}(; parent=ZeroTangent()))
    end

    @testset "to Adjoint" begin
        x = rand(ComplexF64, 3, 4)
        a = adjoint(x)
        ma = collect(a)

        @test a == ProjectTo(a)(ma)
        @test zeros(4, 3) == ProjectTo(a)(ZeroTangent())
        @test zeros(4, 3) == ProjectTo(a)(Tangent{Adjoint}(; parent=ZeroTangent()))
    end

    @testset "to PermutedDimsArray" begin
        a = zeros(3, 5, 4)
        b = PermutedDimsArray(a, (2, 1, 3))
        bc = collect(b)

        @test b == ProjectTo(b)(bc)
        @test b == ProjectTo(b)(ZeroTangent())
    end

    @testset "to SubArray" begin
        x = rand(3, 4)
        sa = view(x, :, 1:2)
        m = collect(sa)

        # make sure it converts the view to the parent type
        @test ProjectTo(sa)(m) isa Matrix
        @test zeros(3, 2) == ProjectTo(sa)(ZeroTangent())
        @test ProjectTo(sa)(Tangent{SubArray}(; parent=ZeroTangent())) isa Matrix
    end
end


text:
"
    Returns spin operators Sx,Sy,Sz,Id for spin s
"
function spinmatrices(s::Union{Rational{Int},Int},elt = ComplexF64)
    N = Int(2*s)

    Sx = zeros(elt,N+1,N+1)
    Sy = zeros(elt,N+1,N+1)
    Sz = zeros(elt,N+1,N+1)

    for row=1:(N+1)
        for col=1:(N+1)
            term=sqrt((s+1)*(row+col-1)-row*col)/2.0

            if (row+1==col)
                Sx[row,col]+=term
                Sy[row,col]-=1im*term
            end

            if(row==col+1)
                Sx[row,col]+=term
                Sy[row,col]+=1im*term
            end

            if(row==col)
                Sz[row,col]+=s+1-row
            end

        end
    end
    return Sx,Sy,Sz,one(Sx)
end

function nonsym_spintensors(s)
    (Sxd,Syd,Szd) = spinmatrices(s)
    sp = ComplexSpace(size(Sxd,1))

    Sx = TensorMap(Sxd,sp,sp);
    Sy = TensorMap(Syd,sp,sp);
    Sz = TensorMap(Szd,sp,sp);

    return Sx,Sy,Sz,one(Sx)
end

"""
bosonic creation anihilation operators with a cutoff
cutoff = maximal number of bosons at one location
"""
function nonsym_bosonictensors(cutoff::Int,elt=ComplexF64)
    creadat = zeros(elt,cutoff+1,cutoff+1);

    for i in 1:cutoff
        creadat[i+1,i] = sqrt(i);
    end

    a⁺ = TensorMap(creadat,ℂ^(cutoff+1),ℂ^(cutoff+1));
    a⁻ = TensorMap(collect(creadat'),ℂ^(cutoff+1),ℂ^(cutoff+1));
    return (a⁺,a⁻)
end

ext: jl

max_stars_repo_name: maartenvd/MPSKitModels

max_issues_repo_name: maartenvd/MPSKitModels

max_forks_repo_name: maartenvd/MPSKitModels


text:
include("mp-datastructures.jl")
include("mp-functions.jl")
using Plots
pyplot()

"""    emptyboard(ngames, nmoves)
Player just loops around an empty board. Returns move number and money as n x 2 array.
Move numbers - out[:,1]
Money - out[:,2]"""
function emptyboard(ngames, nmoves)
    n = ngames * nmoves
    out = zeros(n,2)
    i = 0
    for game in 1:ngames
        mygame = newgame(1500)
        # random starting point
        currentplayer = mygame.players[mygame.currentplayer]
        currentplayer.location = rand(1:40)
        for move in 1:nmoves
            i += 1
            currentplayer = mygame.players[mygame.currentplayer]
            visited = taketurn!(mygame)
            out[i,:] = [move, currentplayer.money]
        end
    end
    return out
end

"""    linearfit(x, y)
Computes least squares fit straight line between scattered x and y points.
Returns pair of points as an x-array and y-array that can be used to directly plot a straight line"""
function linearfit(x, y)
    # sum_xx is sum of squares of x, sum_xy is sum of x*y
    sum_xx = sum(x -> x^2, x)
    sum_xy = sum(x .* y)
    n = length(x)
    # y = bx + a where b is slope and a is intercept
    slope = (n*sum_xy - sum(x)*sum(y)) /
            (n*sum_xx - sum(x)^2)
    intercept = (sum(y) - slope*sum(x)) / n
    xpoints = [minimum(x), maximum(x)]
    ypoints = slope .* xpoints .+ intercept
    return xpoints, ypoints
end

@time balances = emptyboard(200, 80)
x = balances[:,1]
y = balances[:,2]
xtrend, ytrend = linearfit(x, y)
plot(x, y
    ,seriestype=:scatter, seriesalpha=.08, seriescolor="brown"
    ,markerstrokewidth=0, markersize=4.5
    ,gridalpha=.5, legend=false
    ,title="Just Going Around The Empty Board"
    ,xlabel="Move Number", ylabel="Balance (£)"
    ,tickfontsize=10)
plot!(xtrend, ytrend, seriescolor="blue", linewidth=2)


text:
#!/usr/bin/env -S julia -iq
#
# Jerzy Wroczyński (nr. indeksu 250075)
# Obliczenia Naukowe, Laboratorium
# Lista 4., Zadanie 4. — podstawowy test
#

include("interpolation.jl")
using .Interpolation

Interpolation.rysujNnfx(x->sin(x) + cos(x), 0.0, 5.0, 5)

println("end the interactive session to close all plots")

max_stars_repo_name: jerry-sky/academic-notebook

max_issues_repo_name: jerry-sky/academic-notebook

max_forks_repo_name: jerry-sky/academic-notebook


text:
# run_all for unconstrained
function run_all(obj, x0, approx, globals)
    res = []
    _res = minimize!(obj, copy(x0), approx)
    push!(res, _res)
    for glob in globals
      _res = minimize!(obj, copy(x0), (approx, glob))
      printstyled(_res; color=maximum(abs, _res[3]) ≤ sqrt(eps(Float64)) ? :green : :red);println()
      push!(res, _res)
    end
    res
end


text:
# This file is a part of JuliaFEM.
# License is MIT: see https://github.com/JuliaFEM/FEMBasis.jl/blob/master/LICENSE

mutable struct NSolid <: AbstractBasis
    order_u :: Int
    order_v :: Int
    order_w :: Int
    knots_u :: Vector{Float64}
    knots_v :: Vector{Float64}
    knots_w :: Vector{Float64}
    weights :: Array{Float64, 3}
end

function NSolid()
    NSolid(1, 1, 1,
        [-1.0, -1.0, 1.0, 1.0],
        [-1.0, -1.0, 1.0, 1.0],
        [-1.0, -1.0, 1.0, 1.0],
        ones(2, 2, 2))
end

function length(basis::NSolid)
    nu = length(basis.knots_u) - basis.order_u - 1
    nv = length(basis.knots_v) - basis.order_v - 1
    nw = length(basis.knots_w) - basis.order_w - 1
    return nu*nv*nw
end

function size(basis::NSolid)
    return (3, length(basis))
end

function eval_basis!(basis::NSolid, N::Matrix, xi)
    pu = basis.order_u
    pv = basis.order_v
    pw = basis.order_w
    tu = basis.knots_u
    tv = basis.knots_v
    tw = basis.knots_w
    weights = basis.weights
    nu = length(tu)-pu-1
    nv = length(tv)-pv-1
    nw = length(tw)-pw-1
    u, v, w = xi
    n = 1
    for i=1:nu
        for j=1:nv
            for k=1:nw
                A = NURBS(i,pu,u,tu)
                B = NURBS(j,pv,v,tv)
                C = NURBS(k,pw,w,tw)
                N[n] = weights[i,j,k]*A*B*C
                n += 1
            end
        end
    end
    N ./= sum(N)
    return N
end


text:
using Statistics
using Plots

rectangle(xmin, xmax, ymin, ymax) = Plots.Shape([xmin,xmax,xmax,xmin], [ymin,ymin,ymax,ymax])
function Plots_boxbars!(p, x, dx, Q_; boxColor=:lightgray)
    #verrorbars(x, Q_[:,3], Q_[:,1], Q_[:,5])
    for k in 1:length(x)
        x_ = x[k]
        Plots.plot!(p, [x_,x_], Q_[k,[1,5]], color="black")
        Plots.plot!(p, rectangle(x_-dx,x_+dx,Q_[k,2],Q_[k,4]), color=boxColor)
        Plots.plot!(p, [x_-dx,x_+dx], Q_[k,[3,3]], color="black")
        Plots.plot!(p, [x_-dx,x_+dx], Q_[k,[1,1]], color="black")
        Plots.plot!(p, [x_-dx,x_+dx], Q_[k,[5,5]], color="black")
    end
    nothing
end

function Plots_all_estimates!(p, xs, E; group_width=0.3, box_width=0.8,
    max_y=Inf, only_mean=false, markershape=:xcross, showmean=true,
    box_colors = (:lightgray, :lightyellow, :lightgreen, :pink, :lightblue, :cyan))
    n_colors = length(box_colors)
    function limits_tune(xmin, xmax; gap=0.05)
        local dx = xmax - xmin
        (xmin - dx*gap, xmax + dx*gap)
    end

    n_xs = length(xs)
    n_ns = size(E)[2]
    dx = (xs[2]-xs[1])*group_width/n_xs
    xmin = xs[1]-dx*n_ns/2; xmax = xs[end]+dx*(n_ns+2)/2
    ymin = Inf; ymax = -Inf
    # For now: just plug in max cost
    #E[E.==Inf] .= max_y
    Q_ = Vector(undef, n_ns)
    M_ = Vector(undef, n_ns)
    for k in 1:n_ns
        M_[k] = mapslices(x -> mean(x), E[:,k,:], dims=2)
        if ~only_mean
            Q_[k] = mapslices(x -> quantile(x, [0.05,0.25,0.5,0.75,0.95]), E[:,k,:], dims=2)
            ymin = min(ymin, minimum(Q_[k]))
            ymax = max(ymax, maximum(Q_[k]))
        end
    end
    ymax = min(max_y, ymax)

    Plots.plot!(p, xlim=(xmin-n_xs/2*dx, xmax), ylim=(0, ymax), #ylim=limits_tune(ymin, ymax),
    framestyle=:axes, gridlinewidht=3.0, gridstyle=:solid,
    gridalpha=0.2, xticks=xs)
    display(p)
    #result = input("foo")
    if ~only_mean
        for k in 1:n_xs
            xs_ = xs[k]
            Plots.vspan!(p, [xs_-dx*n_xs*1.3,xs_], color=:gray, fillalpha=0.3, linealpha=0)
        end
        for k in 1:n_ns
            col = box_colors[rem(k-1,n_colors)+1]
            xs_ = xs .+ (k .- .5 .- n_ns/2)*dx
            Plots_boxbars!(p, xs_, dx/2*box_width, Q_[k]; boxColor=col)
        end
    end
    if showmean
    for k in 1:n_ns
        xs_ = xs .+ (k .- .5 .- n_ns/2)*dx
        Plots.plot!(p, xs_, M_[k]; seriestype=:scatter, markershape = markershape,
        markercolor = :black, markersize=3)
    end
    end
    #for k in 1:n_xs
    #    col = box_colors[rem(k-1,n_colors)+1]
    #    xs_ = xs[k] .+ (collect(1:n_ns) .- n_ns/2)*dx
    #    Plots_boxbars!(p, xs_, dx/2*box_width, Q_[k]; boxColor=col)
        #M_ = vec(mapslices(x->mean(x), E[k,:,:], dims=2))
        #polymarker(xs_, M_)
    #end
    Plots.plot!(p, legend=false)
    nothing
end

function add_legend!(p, Ns, x_, dx, y, box_colors)
    n_Ns = length(Ns)
    dy = diff(y)[1]/n_Ns
    for i in 1:n_Ns
        y_ = y[2] - dy*(i-.5)
        Plots.plot!(p, rectangle(x_,x_+dx,y_-dy*.4, y_+dy*.4), color=box_colors[i])
        Plots.annotate!(p, x_+1.5*dx, y_, text(string("N=", Ns[i]), :left, 7))
    end
end

function plotCost(ts, AS_iter, BS_iter, Ns, maxIter, max_y=75)
    # Calculate 'normalised cost' = iter*N_particles/length_of_series
    n_ts = length(ts); n_Ns = length(Ns); n_rep = size(AS_iter)[3]
    BS_cost = zeros(n_ts, n_Ns, n_rep)
    AS_cost = zeros(n_ts, n_Ns, n_rep)
    BS_capped = zeros(n_ts, n_Ns)
    AS_capped = zeros(n_ts, n_Ns)
    BS_capval = zeros(n_ts, n_Ns)*NaN
    AS_capval = zeros(n_ts, n_Ns)*NaN
    for i in 1:n_ts
      t = ts[i]
      for j in 1:n_Ns
          N = Ns[j]
          maxCost_ = maxIter[i,j]*N/t
          BS_capped_ = isinf.(BS_iter[i,j,:])
          AS_capped_ = isinf.(AS_iter[i,j,:])
          BS_capped[i,j] = sum(BS_capped_)
          AS_capped[i,j] = sum(AS_capped_)
          BS_cost[i,j,:] .= BS_iter[i,j,:]*N/t
          AS_cost[i,j,:] .= AS_iter[i,j,:]*N/t
          BS_cost[i,j,BS_capped_] .= maxCost_
          AS_cost[i,j,AS_capped_] .= maxCost_
          if any(BS_capped_)
              BS_capval[i,j] = min(maxCost_, max_y)
          end
          if any(AS_capped_)
              AS_capval[i,j] = min(maxCost_, max_y)
          end
      end
    end
    #BS_cost = mapslices(x -> x./ts, mapslices(x -> x.*Ns, BS_iter, dims=2), dims=1)
    #AS_cost = mapslices(x -> x./ts, mapslices(x -> x.*Ns, AS_iter, dims=2), dims=1)
    #maxCost = mapslices(x -> x./ts, mapslices(x -> x.*Ns, maxIter, dims=2), dims=1)
    box_colors = (:lightgray, :cyan, :lightyellow, :lightgreen, :pink, :lightblue)

    p = Plots.plot()
    Plots_all_estimates!(p, ts, [BS_cost AS_cost], showmean=false, box_colors=box_colors)
    add_legend!(p, Ns, 4900,80, [45,75], box_colors)
    Plots_all_estimates!(p, ts, [BS_capval AS_capval]; only_mean=true, markershape=:star8)
    Plots.plot!(p, ylim=(0,max_y), xlabel="T", ylabel="Normalised cost")
    p, AS_capped, BS_capped
end

using JLD2
@load "simpleDemoCompareCostAll_results.jld2"
p, AS_capped, BS_capped = plotCost(ts, AS_iter, BS_iter, Ns, maxIter)
Plots.plot!(p, size=(700,200))
Plots.savefig("simpleDemoCompareCost.pdf")


text:
##########################################################################
# Copyright 2017 Samuel Ridler.
#
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
#     http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.
##########################################################################

# crop raster borders to fit in map borders, return true if successful, false otherwise
# mutates: raster
function cropRaster!(raster::Raster, map::Map)

	# find x, y, and z values to keep
	ix = findall((map.xMin + raster.dx / 2 .<= raster.x) .* (raster.x .<= map.xMax - raster.dx / 2))
	iy = findall((map.yMin + raster.dy / 2 .<= raster.y) .* (raster.y .<= map.yMax - raster.dy / 2))
	x = raster.x[ix]
	y = raster.y[iy]
	z = raster.z[ix,iy]

	# create cropped raster
	croppedRaster = nothing # init
	try
		croppedRaster = Raster(x, y, z)
	catch
		println("Raster does not overlap enough with map; raster not cropped")
		@show length(x), length(y)
		return false
	end

	# overwrite raster field values in raster with those in croppedRaster
	for fname in fieldnames(Raster)
		setfield!(raster, fname, getfield(croppedRaster, fname))
	end

	return true
end

# given raster and x and y indices (i, j) for a cell,
# generate a random location within the cell
# (x[i] +/- 0.5*dx, y[j] +/- 0.5*dy)
function rasterCellRandLocation(raster::Raster, i::Int, j::Int;
	rng::AbstractRNG = GLOBAL_RNG)
	randLocation = Location()
	randLocation.x = raster.x[i] + (rand(rng) - 0.5) * raster.dx
	randLocation.y = raster.y[j] + (rand(rng) - 0.5) * raster.dy
	return randLocation
end

# function rasterZIndexToXIndex(raster::Raster, i::Int)
	# # return index of raster.x for given index of raster.z
	# return mod(i-1, raster.nx) + 1 # x index
# end

# function rasterZIndexToYIndex(raster::Raster, i::Int)
	# # return index of raster.y for given index of raster.z
	# return cld(i, raster.nx) # y index
# end

function rasterZIndexToXYIndices(raster::Raster, i::Int)
	# return (x,y) indices for raster.z[i]
	@assert(1 <= i && i <= length(raster.z))
	c = CartesianIndices(raster.z)[i]
	return c[1], c[2]
end

# generate n random locations for a given raster sampler
function rasterRandLocations(rasterSampler::RasterSampler, n::Int)
	randLocations = Vector{Location}(undef, n)
	zis = rand(rasterSampler.cellDistrRng, n) # z indices of random locations
	for i = 1:n
		zi = zis[i]
		(xi, yi) = rasterZIndexToXYIndices(rasterSampler.raster, zi)
		randLocations[i] = rasterCellRandLocation(rasterSampler.raster, xi, yi; rng = rasterSampler.cellLocRng)
	end
	return randLocations
end

function printRasterSize(raster::Raster)
	println("nx * ny = ", raster.nx * raster.ny)
	println("xMin; xMax; yMin; yMax:")
	println(minimum(raster.x) - raster.dx/2, "; ", maximum(raster.x) + raster.dx/2, "; ",
		minimum(raster.y) - raster.dy/2, "; ", maximum(raster.y) + raster.dy/2)
end


text:
module YaoClifford

using Yao
using QuantumClifford
# using YaoExtensions
using BitBasis
using LinearAlgebra

export StabilizerReg,
    zero_stabilizer_state,
    rand_stabilizer_state,
    PauliString

include("stabilizer_reg.jl")
include("instruct.jl")
include("measure.jl")

end


text:
"""
Type for inversion variable
* `x` : inversion vector that is optimized
* `m` : model without precon
* `gx` : allocation to store gradient of x
"""
mutable struct X{T,N}
	x::Vector{T}
	last_x::Vector{T}  # last_x needed for optim
	lower_x::Vector{T}
	upper_x::Vector{T}
	m::Array{T,N}
	gx::Vector{T}
	gm::Vector{Array{T,N}}
	prior::Vector{T} # prior on x
	precon::Vector{T}
	preconI::Vector{T} # inv of precon
	w::Vector{T} # only used to compute the distance from prior
end

function X(nx::Int, nobj::Int, m=nothing)
	x=zeros(nx)
	last_x=zero(x)
	lower_x=zero(x)
	upper_x=zero(x)
	gx=zeros(nx)
	prior=zeros(nx)

	precon=ones(nx)
	preconI=ones(nx)

	w=ones(nx)

	if(m===nothing)
		m=zeros(nx)
	end
	gm=[zero(m) for iobj in 1:nobj]

	xx=X(x,last_x, upper_x, lower_x, m,gx,gm,prior,precon,preconI,w)

	update_preconI!(xx)
	return xx
end

"""
update inversion variable using model
Note that `m` can be an array of length same as `x`.
"""
function update_x!(X)
	for i in eachindex(X.x)
		@inbounds X.x[i]=X.m[i]*X.precon[i]
	end
end

function update_m!(X)
	for i in eachindex(X.x)
		@inbounds X.m[i]=X.x[i]*X.preconI[i]
	end
end

"""
update ∂ₓJ using ∂ₘJ
"""
function update_gx!(X, α)
	gx=X.gx
	for i in eachindex(gx)
		gx[i]=zero(eltype(gx))
	end
	for (j,gm) in enumerate(X.gm)
		for i in eachindex(gx)
			@inbounds gx[i]+=gm[i]*X.preconI[i]*α[j]
		end
	end
end

function initialize!(X::X)
	T=eltype(X.x)
	for i in eachindex(X.x)
		X.x[i]=zero(T)
		X.m[i]=zero(T)
		X.gx[i]=zero(T)
		X.prior[i]=zero(T)
		X.precon[i]=one(T)
		X.preconI[i]=one(T)
		X.w[i]=one(T)
	end
	for j in 1:length(X.gm)
		X.gm[j][:]=zero(T)
	end
end

function update_precon!(X, precon)
	p=X.precon
	if(precon===nothing)
		for i in eachindex(p)
			p[i]=one(p)
		end
	else
		copyto!(p, precon)
	end
	update_preconI(X)
end
function update_preconI!(X)
	for i in eachindex(X.precon)
		if(iszero(X.precon[i]))
			X.preconI[i]=X.precon[i]
		else
			X.preconI[i]=inv(X.precon[i])
		end
	end
end

function update_weights!(X, weights)
	w=X.w
	if(weights===nothing)
		for i in eachindex(w)
			w[i]=one(eltype(w))
		end
	else
		copyto!(w, weights)
	end
end


text:
using Pkg

Pkg.activate(".")

using CSV, DataFrames
using GLMakie, Colors, ColorSchemes
using FileIO, Downloads
using GeometryBasics

# Manually downloaded the image
# https://visibleearth.nasa.gov/collection/1484/blue-marble
earth_img = load("data/eo_base_2020_clean_3600x1800.png")

# Converting spherical cordinates to x, y, z

function sphere(; r = 1.0, n=32)
    theta = LinRange(0, pi, n)
    phi = LinRange(-pi, pi, 2*n)
    x = [r * cos(phi) * sin(theta) for theta in theta, phi in phi]
    y = [r * sin(phi) * sin(theta) for theta in theta, phi in phi]
    z = [r * cos(theta) for theta in theta, phi in phi]
    return (x, y, z)
end

# Read the site data
# Compiled from two different files: Checked on 11/02/2022
# Location file https://fluxnet.org/sites/site-list-and-pages/
# TimePeriod Available https://fluxnet.org/data/data-availability/

site_data = DataFrame(CSV.File("data/site_data.csv"))

# Function to convert lat lon to the x,y,z coordinate

function toCartesian(lon , lat ; r = 1.02, cxyz = (0,0,0))
    x = cxyz[1] + (r + 1500_000)*cosd(lat)*cosd(lon)
    y = cxyz[2] + (r + 1500_000)*cosd(lat)*sind(lon)
    z = cxyz[3] + (r + 1500_000)*sind(lat)

    return (x, y, z) ./ 1500_000

end

# cylinderMesh = Sphere(Point3f(0.0), 1.0)

lons , lats = site_data.LOCATION_LONG, site_data.LOCATION_LAT

# Magnitude of markers depedning on years available
mag =  site_data.END_YEAR - site_data.STRT_YEAR

toPoints3D = [Point3f([toCartesian(lons[i], lats[i]; r = 0)...]) for i in 1:length(lons)]

# rectMesh = Rect3D(Vec3f(-0.5), Vec3f(2))

set_theme!(theme_minimal())

fig =  Figure(resolution = (800,800), fontsize = 24)

ax = LScene(fig[1,1], show_axis = false)

pltobj = meshscatter!(ax, toPoints3D;
        markersize = mag/2000,
        color =mag,
        colormap = to_colormap(:tab10),
        shading = true,
        ambient = Vec3f(1.0,1.0,1.0),
        )

surface!(ax, sphere(; r = 1.0)...,
        color = tuple.(earth_img, 0.9),
        shading = false,
        transparency = true
        )

Colorbar(fig[1,2], pltobj, label="Number of Years", height = Relative(1/2))
Label(fig[1, 1, Bottom()], "Visualization by @mohitanand")
Label(fig[1, 1, Top()], "Fluxnet data availablity : 211 sites")

zoom!(ax.scene, cameracontrols(ax.scene), 0.65)
rotate!(ax.scene, 3.0)
    ## display(fig)
record(fig, joinpath(@__DIR__, "cflux_white_fast_small.mp4"), framerate = 24) do io
    for i in 3.0:0.015:9.5
        rotate!(ax.scene, i)
        recordframe!(io)  # record a new frame
    end
end
set_theme!()

max_stars_repo_name: melioristic/carb-vis

max_issues_repo_name: CUPofTEAproject/carb-vis

max_forks_repo_name: CUPofTEAproject/carb-vis


text:
# Note: Functions surrounded by a comment blocks are there because `Vararg` is still allocating.
# When Vararg is fast enough, they can simply be removed.

#######################
# UndefBlocksInitializer #
#######################

"""
    UndefBlocksInitializer

Singleton type used in block array initialization, indicating the
array-constructor-caller would like an uninitialized block array. See also
undef_blocks (@ref), an alias for UndefBlocksInitializer().

Examples

≡≡≡≡≡≡≡≡≡≡
```julia
julia> BlockArray(undef_blocks, Matrix{Float32}, [1,2], [3,2])
2×2-blocked 3×5 BlockArray{Float32,2}:
 #undef  #undef  #undef  │  #undef  #undef
 ────────────────────────┼────────────────
 #undef  #undef  #undef  │  #undef  #undef
 #undef  #undef  #undef  │  #undef  #undef
 ```
"""
struct UndefBlocksInitializer end

"""
    undef_blocks

Alias for UndefBlocksInitializer(), which constructs an instance of the singleton
type UndefBlocksInitializer (@ref), used in block array initialization to indicate the
array-constructor-caller would like an uninitialized block array.

Examples

≡≡≡≡≡≡≡≡≡≡
```julia
julia> BlockArray(undef_blocks, Matrix{Float32}, [1,2], [3,2])
2×2-blocked 3×5 BlockArray{Float32,2}:
 #undef  #undef  #undef  │  #undef  #undef
 ------------------------┼----------------
 #undef  #undef  #undef  │  #undef  #undef
 #undef  #undef  #undef  │  #undef  #undef
 ```
"""
const undef_blocks = UndefBlocksInitializer()

##############
# BlockArray #
##############

function _BlockArray end

"""
    BlockArray{T, N, R<:AbstractArray{<:AbstractArray{T,N},N}, BS<:AbstractBlockSizes{N}} <: AbstractBlockArray{T, N}

A `BlockArray` is an array where each block is stored contiguously. This means that insertions and retrieval of blocks
can be very fast and non allocating since no copying of data is needed.

In the type definition, `R` defines the array type that holds the blocks, for example `Matrix{Matrix{Float64}}`.
"""
struct BlockArray{T, N, R <: AbstractArray{<:AbstractArray{T,N},N}, BS<:AbstractBlockSizes{N}} <: AbstractBlockArray{T, N}
    blocks::R
    block_sizes::BS

    global function _BlockArray(blocks::R, block_sizes::BS) where {T, N, R<:AbstractArray{<:AbstractArray{T,N},N}, BS<:AbstractBlockSizes{N}}
        new{T, N, R, BS}(blocks, block_sizes)
    end
end

# Auxilary outer constructors
function _BlockArray(blocks::R, block_sizes::Vararg{AbstractVector{Int}, N}) where {T, N, R<:AbstractArray{<:AbstractArray{T,N},N}}
    return _BlockArray(blocks, BlockSizes(block_sizes...))
end

function _BlockArray(blocks::R, block_sizes::BS) where {T, N, R<:AbstractArray{<:AbstractArray{V,N} where V,N}, BS<:AbstractBlockSizes{N}}
    _BlockArray(convert(AbstractArray{AbstractArray{mapreduce(eltype,promote_type,blocks),N},N}, blocks), block_sizes)
end

const BlockMatrix{T, R <: AbstractMatrix{<:AbstractMatrix{T}}} = BlockArray{T, 2, R}
const BlockVector{T, R <: AbstractVector{<:AbstractVector{T}}} = BlockArray{T, 1, R}
const BlockVecOrMat{T, R} = Union{BlockMatrix{T, R}, BlockVector{T, R}}

################
# Constructors #
################

@inline _BlockArray(::Type{R}, block_sizes::Vararg{AbstractVector{Int}, N}) where {T, N, R<:AbstractArray{<:AbstractArray{T,N},N}} =
    _BlockArray(R, BlockSizes(block_sizes...))

function _BlockArray(::Type{R}, block_sizes::BlockSizes{N}) where {T, N, R<:AbstractArray{<:AbstractArray{T,N},N}}
    n_blocks = nblocks(block_sizes)
    blocks = R(undef, n_blocks)
    _BlockArray(blocks, block_sizes)
end

@inline undef_blocks_BlockArray(::Type{R}, block_sizes::Vararg{AbstractVector{Int}, N}) where {T, N, R<:AbstractArray{<:AbstractArray{T,N},N}} =
    _BlockArray(R, block_sizes...)

"""
Constructs a `BlockArray` with uninitialized blocks from a block type `R` with sizes defind by `block_sizes`.

```jldoctest; setup = quote using BlockArrays end
julia> BlockArray(undef_blocks, Matrix{Float64}, [1,3], [2,2])
2×2-blocked 4×4 BlockArray{Float64,2}:
 #undef  │  #undef  #undef  #undef  │
 --------┼--------------------------┼
 #undef  │  #undef  #undef  #undef  │
 #undef  │  #undef  #undef  #undef  │
 --------┼--------------------------┼
 #undef  │  #undef  #undef  #undef  │
```
"""
@inline BlockArray(::UndefBlocksInitializer, ::Type{R}, block_sizes::Vararg{AbstractVector{Int}, N}) where {T, N, R<:AbstractArray{T,N}} =
    undef_blocks_BlockArray(Array{R,N}, block_sizes...)

@inline BlockArray{T}(::UndefBlocksInitializer, block_sizes::Vararg{AbstractVector{Int}, N}) where {T, N} =
    BlockArray(undef_blocks, Array{T,N}, block_sizes...)

@inline BlockArray{T,N}(::UndefBlocksInitializer, block_sizes::Vararg{AbstractVector{Int}, N}) where {T, N} =
    BlockArray(undef_blocks, Array{T,N}, block_sizes...)

@inline BlockArray{T,N,R}(::UndefBlocksInitializer, block_sizes::Vararg{AbstractVector{Int}, N}) where {T, N, R<:AbstractArray{<:AbstractArray{T,N},N}} =
    undef_blocks_BlockArray(R, block_sizes...)

@generated function initialized_blocks_BlockArray(::Type{R}, block_sizes::BlockSizes{N}) where R<:AbstractArray{V,N} where {T,N,V<:AbstractArray{T,N}}
    return quote
        block_arr = _BlockArray(R, block_sizes)
        @nloops $N i i->(1:nblocks(block_sizes, i)) begin
            block_index = @ntuple $N i
            setblock!(block_arr, similar(V, blocksize(block_sizes, block_index)), block_index...)
        end

        return block_arr
    end
end

initialized_blocks_BlockArray(::Type{R}, block_sizes::Vararg{AbstractVector{Int}, N}) where {T, N, R<:AbstractArray{<:AbstractArray{T,N},N}} =
    initialized_blocks_BlockArray(R, BlockSizes(block_sizes...))

@inline BlockArray{T}(::UndefInitializer, block_sizes::BlockSizes{N}) where {T, N} =
    initialized_blocks_BlockArray(Array{Array{T,N},N}, block_sizes)

@inline BlockArray{T, N}(::UndefInitializer, block_sizes::BlockSizes{N}) where {T, N} =
    initialized_blocks_BlockArray(Array{Array{T,N},N}, block_sizes)

@inline BlockArray{T, N, R}(::UndefInitializer, block_sizes::BlockSizes{N}) where {T, N, R<:AbstractArray{<:AbstractArray{T,N},N}} =
    initialized_blocks_BlockArray(R, block_sizes)

@inline BlockArray{T}(::UndefInitializer, block_sizes::Vararg{AbstractVector{Int}, N}) where {T, N} =
    initialized_blocks_BlockArray(Array{Array{T,N},N}, block_sizes...)

@inline BlockArray{T, N}(::UndefInitializer, block_sizes::Vararg{AbstractVector{Int}, N}) where {T, N} =
    initialized_blocks_BlockArray(Array{Array{T,N},N}, block_sizes...)

@inline BlockArray{T, N, R}(::UndefInitializer, block_sizes::Vararg{AbstractVector{Int}, N}) where {T, N, R<:AbstractArray{<:AbstractArray{T,N},N}} =
    initialized_blocks_BlockArray(R, block_sizes...)

function BlockArray(arr::AbstractArray{T, N}, block_sizes::Vararg{AbstractVector{Int}, N}) where {T,N}
    for i in 1:N
        if sum(block_sizes[i]) != size(arr, i)
            throw(DimensionMismatch("block size for dimension $i: $(block_sizes[i]) does not sum to the array size: $(size(arr, i))"))
        end
    end
    BlockArray(arr, BlockSizes(block_sizes...))
end

@generated function BlockArray(arr::AbstractArray{T, N}, block_sizes::BlockSizes{N}) where {T,N}
    return quote
        block_arr = _BlockArray(Array{typeof(arr),N}, block_sizes)
        @nloops $N i i->(1:nblocks(block_sizes, i)) begin
            block_index = @ntuple $N i
            indices = globalrange(block_sizes, block_index)
            setblock!(block_arr, arr[indices...], block_index...)
        end

        return block_arr
    end
end

BlockVector(blocks::AbstractVector, block_sizes::AbstractBlockSizes{1}) = BlockArray(blocks, block_sizes)
BlockVector(blocks::AbstractVector, block_sizes::AbstractVector{Int}) = BlockArray(blocks, block_sizes)
BlockMatrix(blocks::AbstractMatrix, block_sizes::AbstractBlockSizes{2}) = BlockArray(blocks, block_sizes)
BlockMatrix(blocks::AbstractMatrix, block_sizes::Vararg{AbstractVector{Int},2}) = BlockArray(blocks, block_sizes...)

"""
    mortar(blocks::AbstractArray)
    mortar(blocks::AbstractArray{R, N}, sizes_1, sizes_2, ..., sizes_N)
    mortar(blocks::AbstractArray{R, N}, block_sizes::BlockSizes{N})

Construct a `BlockArray` from `blocks`.  `block_sizes` is computed from
`blocks` if it is not given.

# Examples
```jldoctest; setup = quote using BlockArrays end
julia> blocks = permutedims(reshape([
                  1ones(1, 3), 2ones(1, 2),
                  3ones(2, 3), 4ones(2, 2),
              ], (2, 2)))
2×2 Array{Array{Float64,2},2}:
 [1.0 1.0 1.0]               [2.0 2.0]
 [3.0 3.0 3.0; 3.0 3.0 3.0]  [4.0 4.0; 4.0 4.0]

julia> mortar(blocks)
2×2-blocked 3×5 BlockArray{Float64,2}:
 1.0  1.0  1.0  │  2.0  2.0
 ───────────────┼──────────
 3.0  3.0  3.0  │  4.0  4.0
 3.0  3.0  3.0  │  4.0  4.0

julia> ans == mortar(
                  (1ones(1, 3), 2ones(1, 2)),
                  (3ones(2, 3), 4ones(2, 2)),
              )
true
```
"""
mortar(blocks::AbstractArray{R, N}, block_sizes::AbstractBlockSizes{N}) where {R, N} =
    _BlockArray(blocks, block_sizes)

mortar(blocks::AbstractArray{R, N}, block_sizes::Vararg{AbstractVector{Int}, N}) where {R, N} =
    _BlockArray(blocks, block_sizes...)

mortar(blocks::AbstractArray) = mortar(blocks, sizes_from_blocks(blocks))

sizes_from_blocks(blocks) = sizes_from_blocks(blocks, axes(blocks)) # allow overriding on axes

function sizes_from_blocks(blocks::AbstractArray{<:Any, N}, _) where N
    if length(blocks) == 0
        return zeros.(Int, size(blocks))
    end
    if !all(b -> ndims(b) == N, blocks)
        error("All blocks must have ndims consistent with ndims = $N of `blocks` array.")
    end
    fullsizes = map!(size, Array{NTuple{N,Int}, N}(undef, size(blocks)), blocks)
    block_sizes = ntuple(ndims(blocks)) do i
        [s[i] for s in view(fullsizes, ntuple(j -> j == i ? (:) : 1, ndims(blocks))...)]
    end
    checksizes(fullsizes, block_sizes)
    return BlockSizes(block_sizes...)
end

getsizes(block_sizes, block_index) = getindex.(block_sizes, block_index)

@generated function checksizes(fullsizes::Array{NTuple{N,Int}, N}, block_sizes::NTuple{N,Vector{Int}}) where N
    quote
        @nloops $N i fullsizes begin
            block_index = @ntuple $N i
            if fullsizes[block_index...] != getsizes(block_sizes, block_index)
                error("size(blocks[", strip(repr(block_index), ['(', ')']),
                      "]) (= ", fullsizes[block_index...],
                      ") is incompatible with expected size: ",
                      getsizes(block_sizes, block_index))
            end
        end
        return fullsizes
    end
end

"""
    mortar((block_11, ..., block_1m), ... (block_n1, ..., block_nm))

Construct a `BlockMatrix` with `n * m`  blocks.  Each `block_ij` must be an
`AbstractMatrix`.
"""
mortar(rows::Vararg{NTuple{M, AbstractMatrix}}) where M =
    mortar(permutedims(reshape(
        foldl(append!, rows, init=eltype(eltype(rows))[]),
        M, length(rows))))

# Convert AbstractArrays that conform to block array interface
convert(::Type{BlockArray{T,N,R}}, A::BlockArray{T,N,R}) where {T,N,R} = A
convert(::Type{BlockArray{T,N}}, A::BlockArray{T,N}) where {T,N} = A
convert(::Type{BlockArray{T}}, A::BlockArray{T}) where {T} = A
convert(::Type{BlockArray}, A::BlockArray) = A

BlockArray{T, N}(A::AbstractArray{T2, N}) where {T,T2,N} =
    BlockArray(Array{T, N}(A), blocksizes(A))
BlockArray{T1}(A::AbstractArray{T2, N}) where {T1,T2,N} = BlockArray{T1, N}(A)
BlockArray(A::AbstractArray{T, N}) where {T,N} = BlockArray{T, N}(A)

convert(::Type{BlockArray{T, N}}, A::AbstractArray{T2, N}) where {T,T2,N} =
    BlockArray(convert(Array{T, N}, A), blocksizes(A))
convert(::Type{BlockArray{T1}}, A::AbstractArray{T2, N}) where {T1,T2,N} =
    convert(BlockArray{T1, N}, A)
convert(::Type{BlockArray}, A::AbstractArray{T, N}) where {T,N} =
    convert(BlockArray{T, N}, A)

copy(A::BlockArray) = _BlockArray(copy.(A.blocks), copy(A.block_sizes))

################################
# AbstractBlockArray Interface #
################################
@inline blocksizes(block_array::BlockArray) = block_array.block_sizes

@inline function getblock(block_arr::BlockArray{T,N}, block::Vararg{Integer, N}) where {T,N}
    @boundscheck blockcheckbounds(block_arr, block...)
    block_arr.blocks[block...]
end

@inline function Base.getindex(block_arr::BlockArray{T,N}, blockindex::BlockIndex{N}) where {T,N}
    @boundscheck blockcheckbounds(block_arr, Block(blockindex.I))
    @inbounds block = getblock(block_arr, blockindex.I...)
    @boundscheck checkbounds(block, blockindex.α...)
    @inbounds v = block[blockindex.α...]
    return v
end

###########################
# AbstractArray Interface #
###########################

@inline function Base.similar(block_array::BlockArray{T,N}, ::Type{T2}) where {T,N,T2}
    _BlockArray(similar(block_array.blocks, Array{T2, N}), copy(blocksizes(block_array)))
end

@inline function Base.getindex(block_arr::BlockArray{T, N}, i::Vararg{Integer, N}) where {T,N}
    @boundscheck checkbounds(block_arr, i...)
    @inbounds v = block_arr[global2blockindex(blocksizes(block_arr), i)]
    return v
end

@inline function Base.setindex!(block_arr::BlockArray{T, N}, v, i::Vararg{Integer, N}) where {T,N}
    @boundscheck checkbounds(block_arr, i...)
    @inbounds block_arr[global2blockindex(blocksizes(block_arr), i)] = v
    return block_arr
end

############
# Indexing #
############

function _check_setblock!(block_arr::BlockArray{T, N}, v, block::NTuple{N, Integer}) where {T,N}
    for i in 1:N
        if size(v, i) != blocksize(block_arr, i, block[i])
            throw(DimensionMismatch(string("tried to assign $(size(v)) array to ", blocksize(block_arr, block), " block")))
        end
    end
end

@inline function setblock!(block_arr::BlockArray{T, N}, v, block::Vararg{Integer, N}) where {T,N}
    @boundscheck blockcheckbounds(block_arr, block...)
    @boundscheck _check_setblock!(block_arr, v, block)
    @inbounds block_arr.blocks[block...] = v
    return block_arr
end

@propagate_inbounds function Base.setindex!(block_array::BlockArray{T, N}, v, block_index::BlockIndex{N}) where {T,N}
    getblock(block_array, block_index.I...)[block_index.α...] = v
end

Base.dataids(arr::BlockArray) = (dataids(arr.blocks)..., dataids(arr.block_sizes)...)
# This is not entirely valid.  In principle, we have to concatenate
# all dataids of all blocks.  However, it makes `dataids` non-inferable.

########
# Misc #
########

@generated function Base.Array(block_array::BlockArray{T, N, R}) where {T,N,R}
    # TODO: This will fail for empty block array
    return quote
        block_sizes = blocksizes(block_array)
        arr = similar(block_array.blocks[1], size(block_array)...)
        @nloops $N i i->(1:nblocks(block_sizes, i)) begin
            block_index = @ntuple $N i
            indices = globalrange(block_sizes, block_index)
            arr[indices...] = getblock(block_array, block_index...)
        end

        return arr
    end
end

@generated function copyto!(block_array::BlockArray{T, N, R}, arr::R) where {T,N,R <: AbstractArray}
    return quote
        block_sizes = blocksizes(block_array)

        @nloops $N i i->(1:nblocks(block_sizes, i)) begin
            block_index = @ntuple $N i
            indices = globalrange(block_sizes, block_index)
            copyto!(getblock(block_array, block_index...), arr[indices...])
        end

        return block_array
    end
end

function Base.fill!(block_array::BlockArray, v)
    for block in block_array.blocks
        fill!(block, v)
    end
    block_array
end

function lmul!(α::Number, block_array::BlockArray)
    for block in block_array.blocks
        lmul!(α, block)
    end
    block_array
end

function rmul!(block_array::BlockArray, α::Number)
    for block in block_array.blocks
        rmul!(block, α)
    end
    block_array
end


text:
module BED

import Automa
import Automa.RegExp: @re_str
import BGZFStreams
import BioCore
import BufferedStreams
import ColorTypes
import FixedPointNumbers: N0f8
import GenomicFeatures

include("record.jl")
include("reader.jl")
include("writer.jl")

end


text:
module TestAqua

using Aqua
using BangBang
using Setfield
using Test

Aqua.test_all(
    BangBang;
    ambiguities=(exclude=[Base.get, Setfield.set, Setfield.modify],),
)

@testset "Compare test/Project.toml and test/environments/main/Project.toml" begin
    @test Text(read(joinpath(@__DIR__, "Project.toml"), String)) ==
          Text(read(joinpath(@__DIR__, "environments", "main", "Project.toml"), String))
end

end  # module


text:
module robustml

#export ftest
#export ftest2
#include("test_func.jl")

export abs_val
export obj
export median_index
export buckets
export cvfolds
export block_X
export block_y
export mom_obj
export f_grad
export block
export block_cv
export soft_threshold
export pg_als
export pg_mom_als
export pg_mom_als_w

include("MOM_ALS.jl")

end

max_stars_repo_name: jstriaukas/robustml

max_issues_repo_name: jstriaukas/robustml

max_forks_repo_name: jstriaukas/robustml


text:
# This file is auto-generated by AWSMetadata.jl
using AWS
using AWS.AWSServices: pi
using AWS.Compat
using AWS.UUIDs

"""
    describe_dimension_keys(end_time, group_by, identifier, metric, service_type, start_time)
    describe_dimension_keys(end_time, group_by, identifier, metric, service_type, start_time, params::Dict{String,<:Any})

For a specific time period, retrieve the top N dimension keys for a metric.  Each response
element returns a maximum of 500 bytes. For larger elements, such as SQL statements, only
the first 500 bytes are returned.

# Arguments
- `end_time`: The date and time specifying the end of the requested time series data. The
  value specified is exclusive, which means that data points less than (but not equal to)
  EndTime are returned. The value for EndTime must be later than the value for StartTime.
- `group_by`: A specification for how to aggregate the data points from a query result. You
  must specify a valid dimension group. Performance Insights returns all dimensions within
  this group, unless you provide the names of specific dimensions within this group. You can
  also request that Performance Insights return a limited number of values for a dimension.
- `identifier`: An immutable, AWS Region-unique identifier for a data source. Performance
  Insights gathers metrics from this data source. To use an Amazon RDS instance as a data
  source, you specify its DbiResourceId value. For example, specify
  db-FAIHNTYBKTGAUSUZQYPDS2GW4A
- `metric`: The name of a Performance Insights metric to be measured. Valid values for
  Metric are:    db.load.avg - a scaled representation of the number of active sessions for
  the database engine.    db.sampledload.avg - the raw number of active sessions for the
  database engine.   If the number of active sessions is less than an internal Performance
  Insights threshold, db.load.avg and db.sampledload.avg are the same value. If the number of
  active sessions is greater than the internal threshold, Performance Insights samples the
  active sessions, with db.load.avg showing the scaled values, db.sampledload.avg showing the
  raw values, and db.sampledload.avg less than db.load.avg. For most use cases, you can query
  db.load.avg only.
- `service_type`: The AWS service for which Performance Insights will return metrics. The
  only valid value for ServiceType is RDS.
- `start_time`: The date and time specifying the beginning of the requested time series
  data. You must specify a StartTime within the past 7 days. The value specified is
  inclusive, which means that data points equal to or greater than StartTime are returned.
  The value for StartTime must be earlier than the value for EndTime.

# Optional Parameters
Optional parameters can be passed as a `params::Dict{String,<:Any}`. Valid keys are:
- `"Filter"`: One or more filters to apply in the request. Restrictions:   Any number of
  filters by the same dimension, as specified in the GroupBy or Partition parameters.   A
  single filter for any other dimension in this dimension group.
- `"MaxResults"`: The maximum number of items to return in the response. If more items
  exist than the specified MaxRecords value, a pagination token is included in the response
  so that the remaining results can be retrieved.
- `"NextToken"`: An optional pagination token provided by a previous request. If this
  parameter is specified, the response includes only records beyond the token, up to the
  value specified by MaxRecords.
- `"PartitionBy"`: For each dimension specified in GroupBy, specify a secondary dimension
  to further subdivide the partition keys in the response.
- `"PeriodInSeconds"`: The granularity, in seconds, of the data points returned from
  Performance Insights. A period can be as short as one second, or as long as one day (86400
  seconds). Valid values are:    1 (one second)    60 (one minute)    300 (five minutes)
  3600 (one hour)    86400 (twenty-four hours)   If you don't specify PeriodInSeconds, then
  Performance Insights chooses a value for you, with a goal of returning roughly 100-200 data
  points in the response.
"""
describe_dimension_keys(EndTime, GroupBy, Identifier, Metric, ServiceType, StartTime; aws_config::AbstractAWSConfig=global_aws_config()) = pi("DescribeDimensionKeys", Dict{String, Any}("EndTime"=>EndTime, "GroupBy"=>GroupBy, "Identifier"=>Identifier, "Metric"=>Metric, "ServiceType"=>ServiceType, "StartTime"=>StartTime); aws_config=aws_config)
describe_dimension_keys(EndTime, GroupBy, Identifier, Metric, ServiceType, StartTime, params::AbstractDict{String}; aws_config::AbstractAWSConfig=global_aws_config()) = pi("DescribeDimensionKeys", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("EndTime"=>EndTime, "GroupBy"=>GroupBy, "Identifier"=>Identifier, "Metric"=>Metric, "ServiceType"=>ServiceType, "StartTime"=>StartTime), params)); aws_config=aws_config)

"""
    get_dimension_key_details(group, group_identifier, identifier, service_type)
    get_dimension_key_details(group, group_identifier, identifier, service_type, params::Dict{String,<:Any})

Get the attributes of the specified dimension group for a DB instance or data source. For
example, if you specify a SQL ID, GetDimensionKeyDetails retrieves the full text of the
dimension db.sql.statement associated with this ID. This operation is useful because
GetResourceMetrics and DescribeDimensionKeys don't support retrieval of large SQL statement
text.

# Arguments
- `group`: The name of the dimension group. The only valid value is db.sql. Performance
  Insights searches the specified group for the dimension group ID.
- `group_identifier`: The ID of the dimension group from which to retrieve dimension
  details. For dimension group db.sql, the group ID is db.sql.id.
- `identifier`: The ID for a data source from which to gather dimension data. This ID must
  be immutable and unique within an AWS Region. When a DB instance is the data source,
  specify its DbiResourceId value. For example, specify db-ABCDEFGHIJKLMNOPQRSTU1VW2X.
- `service_type`: The AWS service for which Performance Insights returns data. The only
  valid value is RDS.

# Optional Parameters
Optional parameters can be passed as a `params::Dict{String,<:Any}`. Valid keys are:
- `"RequestedDimensions"`: A list of dimensions to retrieve the detail data for within the
  given dimension group. For the dimension group db.sql, specify either the full dimension
  name db.sql.statement or the short dimension name statement. If you don't specify this
  parameter, Performance Insights returns all dimension data within the specified dimension
  group.
"""
get_dimension_key_details(Group, GroupIdentifier, Identifier, ServiceType; aws_config::AbstractAWSConfig=global_aws_config()) = pi("GetDimensionKeyDetails", Dict{String, Any}("Group"=>Group, "GroupIdentifier"=>GroupIdentifier, "Identifier"=>Identifier, "ServiceType"=>ServiceType); aws_config=aws_config)
get_dimension_key_details(Group, GroupIdentifier, Identifier, ServiceType, params::AbstractDict{String}; aws_config::AbstractAWSConfig=global_aws_config()) = pi("GetDimensionKeyDetails", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("Group"=>Group, "GroupIdentifier"=>GroupIdentifier, "Identifier"=>Identifier, "ServiceType"=>ServiceType), params)); aws_config=aws_config)

"""
    get_resource_metrics(end_time, identifier, metric_queries, service_type, start_time)
    get_resource_metrics(end_time, identifier, metric_queries, service_type, start_time, params::Dict{String,<:Any})

Retrieve Performance Insights metrics for a set of data sources, over a time period. You
can provide specific dimension groups and dimensions, and provide aggregation and filtering
criteria for each group.  Each response element returns a maximum of 500 bytes. For larger
elements, such as SQL statements, only the first 500 bytes are returned.

# Arguments
- `end_time`: The date and time specifying the end of the requested time series data. The
  value specified is exclusive - data points less than (but not equal to) EndTime will be
  returned. The value for EndTime must be later than the value for StartTime.
- `identifier`: An immutable, AWS Region-unique identifier for a data source. Performance
  Insights gathers metrics from this data source. To use a DB instance as a data source,
  specify its DbiResourceId value. For example, specify db-FAIHNTYBKTGAUSUZQYPDS2GW4A.
- `metric_queries`: An array of one or more queries to perform. Each query must specify a
  Performance Insights metric, and can optionally specify aggregation and filtering criteria.
- `service_type`: The AWS service for which Performance Insights returns metrics. The only
  valid value for ServiceType is RDS.
- `start_time`: The date and time specifying the beginning of the requested time series
  data. You can't specify a StartTime that's earlier than 7 days ago. The value specified is
  inclusive - data points equal to or greater than StartTime will be returned. The value for
  StartTime must be earlier than the value for EndTime.

# Optional Parameters
Optional parameters can be passed as a `params::Dict{String,<:Any}`. Valid keys are:
- `"MaxResults"`: The maximum number of items to return in the response. If more items
  exist than the specified MaxRecords value, a pagination token is included in the response
  so that the remaining results can be retrieved.
- `"NextToken"`: An optional pagination token provided by a previous request. If this
  parameter is specified, the response includes only records beyond the token, up to the
  value specified by MaxRecords.
- `"PeriodInSeconds"`: The granularity, in seconds, of the data points returned from
  Performance Insights. A period can be as short as one second, or as long as one day (86400
  seconds). Valid values are:    1 (one second)    60 (one minute)    300 (five minutes)
  3600 (one hour)    86400 (twenty-four hours)   If you don't specify PeriodInSeconds, then
  Performance Insights will choose a value for you, with a goal of returning roughly 100-200
  data points in the response.
"""
get_resource_metrics(EndTime, Identifier, MetricQueries, ServiceType, StartTime; aws_config::AbstractAWSConfig=global_aws_config()) = pi("GetResourceMetrics", Dict{String, Any}("EndTime"=>EndTime, "Identifier"=>Identifier, "MetricQueries"=>MetricQueries, "ServiceType"=>ServiceType, "StartTime"=>StartTime); aws_config=aws_config)
get_resource_metrics(EndTime, Identifier, MetricQueries, ServiceType, StartTime, params::AbstractDict{String}; aws_config::AbstractAWSConfig=global_aws_config()) = pi("GetResourceMetrics", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("EndTime"=>EndTime, "Identifier"=>Identifier, "MetricQueries"=>MetricQueries, "ServiceType"=>ServiceType, "StartTime"=>StartTime), params)); aws_config=aws_config)


text:
using Serialization, Pkg, DocSeeker

db_path(pkg, env) = joinpath(@__DIR__, "..", "db", string(pkg, "-", hash(env), ".db"))

function create_db(pkg, env)
    sympkg = Symbol(pkg)
    db = db_path(pkg, env)
    cd(env) do
        open(db, "w+") do io
            if pkg == "Base"
                serialize(io, DocSeeker.alldocs(Base))
            else
                mod = Main.eval(quote
                    using $sympkg
                    $sympkg
                end)
                serialize(io, DocSeeker.alldocs(mod))
            end
        end
    end
end

pkg, env = ARGS

Pkg.activate(env)
create_db(pkg, env)


text:
module HessMain
using Hecke
using ..HessQRModule, ..GenericRound2
import ..HessQRModule: HessQR

function GenericRound2.integral_closure(S::HessQR, F::Generic.FunctionField{T}) where {T}
  return GenericRound2._integral_closure(S, F)
end

_lcm(a::fmpq, b::fmpq) = fmpq(lcm(numerator(a), numerator(b)), gcd(denominator(a), denominator(b)))
_gcd(a::fmpq, b::fmpq) = fmpq(gcd(numerator(a), numerator(b)), lcm(denominator(a), denominator(b)))

#=
Hallo Claus,

für R Hauptidealring geht es dann wohl so:

Sei M die Trafomatrix zwischen einer R<t>- und einer K[t]-Basis, wobei
die Zeilen zu R<t> und die Spalten zu K[t] korrespondieren.  Es ist
det(M) neq 0. Die Zeilen von M können durch Multiplikation mit primitiven
Polynomen aus R[t] so skaliert werden, dass M nur Polynome aus R[t] als
Einträge hat. Dies ist eine R<t>-unimodulare Operation.

Wenn det(M) Inhalt eins hat, dann ist M R<t>-unimodular, durch
Basiswechsel beim R<t>-Modul läßt sich also M = Einheitsmatrix
erreichen. Also korrespondieren die Spalten von M zu einer Basis des
Schnitts des K[t]- und R<t>-Moduls.

Wenn det(M) nicht Inhalt eins hat, dann gibt es einen Primteiler p des
Inhalts. Dann ist det(M) mod p = 0, also hat die
Spaltenhermitenormalform von M mod p Nullspalten.  Diese
Spaltenhermitenormalform wird durch durch Multiplikation von rechts mit
Elementarmatrizen über R/p[t] erreicht. Diese Elementarmatrizen können
modulo p geliftet werden, so dass sie Einträge aus R[t] haben und
R[t]-unimodular sind.  Multiplikation von rechts mit diesen gelifteten
Elementarmatrizen liefert eine Matrix, die mod p die oben erwähnte
Spaltenhermitenormalform ist. Insbesondere ist mindestens eine der
hinteren Spalten über R[t] durch p teilbar. Also diese Spalten solange
durch p dividieren, bis keine Spalte mehr durch p dividierbar ist.  Das
ist eine k[t]-unimodulare Operation. Nun ist auch der p-Anteil des
Inhalts von det(M) echt verringert worden. Iteration mit diesem oder
anderen p jenachdem, was noch im Inhalt übrig ist, führt schließlich zu
Inhalt eins.

Viele Grüße,
Florian
=#

function florian(M::MatElem{<:Generic.Rat{fmpq}}, R::FmpqPolyRing, S::HessQR)
  Qt = base_ring(M)
  n = nrows(M)
  #step 1: make integral

  MM, d = integral_split(M, R)
  #M and d are in /over Q[x] (of type Q[x])
  cM = reduce(_gcd, map(content, MM), init = fmpq(1,1))
  MM *= inv(cM) #should now be in Z[x]!
  cd = content(d)
  d *= inv(cd)  #should be in Z[x] and primitive, so

  T1 = identity_matrix(Qt, n)*d
  T2 = identity_matrix(Qt, n)*inv(cM)*inv(cd)

  @assert T1*M*T2 == MM
  @assert isone(integral_split(T1, S)[2])
  @assert isone(integral_split(T2, R)[2])
  @assert isunit(integral_split(det(T1), S)[1])
  @assert isunit(integral_split(det(T2), R)[1])

  de = det(MM)
  for p = keys(factor(content(de), ZZ).fac)
    #step 2: do a HNF mod p and record the lifted operations
    k = GF(p)
    kt = PolynomialRing(k, cached = false)[1]
    while true
      H = map(kt, MM)
      piv = 1
      for i=1:n
        if iszero(H[i,piv])
          j = piv+1
          while j<= n && iszero(H[i, j])
            j += 1
          end
          if j>n
            continue
          end
          H[:, piv], H[:, j] = H[:, j], H[:, piv]
          T2[:, piv], T2[:, j] = T2[:, j], T2[:, piv]
          MM[:, piv], MM[:, j] = MM[:, j], MM[:, piv]
        end
        for j=piv+1:n
          while !iszero(H[i, j])
            q, r = divrem(H[i, j], H[i,piv])
            H[:, j] = H[:, j] - q*H[:, piv]
            @assert H[i, j] == r
            T2[:, j] = T2[:, j] - Qt(Hecke.lift(Hecke.Globals.Zx, q))*T2[:, piv]
            MM[:, j] = MM[:, j] - R(Hecke.lift(Hecke.Globals.Zx, q))*MM[:, piv]
            if iszero(r)
              break
            end
            H[:, piv], H[:, j] = H[:, j], H[:, piv]
            T2[:, piv], T2[:, j] = T2[:, j], T2[:, piv]
            MM[:, piv], MM[:, j] = MM[:, j], MM[:, piv]
          end
        end
        @assert !iszero(H[i,piv])
        #= Don't reduce upwards...it does not change the "being zero or not"
        for j=1:piv-1
          q, r = divrem(H[i, j], H[i,piv])
          H[:, j] = H[:, j] - q*H[:, piv]
          @assert H[i, j] == r
          T2[:, j] = T2[:, j] - Qt(Hecke.lift(Hecke.Globals.Zx, q))*T2[:, i]
          MM[:, j] = MM[:, j] - R(Hecke.lift(Hecke.Globals.Zx, q))*MM[:, i]
        end
        =#
        piv += 1
      end
      done = true
      for i=1:n
        if iszero(H[:,i])
          done = false
          T2[:, i] *= Qt(fmpq(1, p))
          MM[:, i] *= R(fmpq(1, p))
        end
      end
      @assert T1*M*T2 == MM
      done && break
    end
  end
  return M, T1, T2
end

function GenericRound2.integral_closure(Zx::FmpzPolyRing, F::Generic.FunctionField)
  Qt = base_ring(F)
  t = gen(Qt)
  S = HessQR(Zx, Qt)
  R = parent(numerator(t))
  o1 = integral_closure(S, F)
  o2 = integral_closure(R, F)
  if isdefined(o1, :trans)
    T = o1.trans
  else
    T = identity_matrix(Qt, degree(F))
  end
  if isdefined(o2, :itrans)
    T = T * o2.itrans
  end
  _, T1, T2 = florian(T, R, S)

  o3 = GenericRound2.Order(Zx, F, true)
  if isdefined(o2, :trans)
    oo2 = GenericRound2.Order(o3, integral_split(inv(T2)*o2.trans, Zx)..., check = false)
  else
    oo2 = GenericRound2.Order(o3, integral_split(inv(T2), Zx)..., check = false)
  end
  return oo2

  #for testing....
  H, TT1 = hnf_with_transform(map_entries(S, T1*T*T2))
  @assert isone(H)
  T1 = map_entries(Qt, TT1)*T1
  if isdefined(o1, :trans)
    oo1 = GenericRound2.Order(o3, integral_split(T1*o1.trans, Zx)..., check = false)
  else
    oo1 = GenericRound2.Order(o3, integral_split(T1, Zx)..., check = false)
  end
  return oo1, oo2
end

function Base.denominator(a::Generic.Rat{fmpq}, S::FmpzPolyRing)
  return integral_split(a, S)[2]
end

function Base.numerator(a::Generic.Rat{fmpq}, S::FmpzPolyRing)
  return integral_split(a, S)[1]
end

function Hecke.integral_split(a::Generic.Rat{fmpq}, S::FmpzPolyRing)
  #TODO: feels too complicated....
  if iszero(a)
    return zero(S), one(S)
  end
  n = numerator(a)
  d = denominator(a)
  dn = reduce(lcm, map(denominator, coefficients(n)), init = fmpz(1))
  dd = reduce(lcm, map(denominator, coefficients(d)), init = fmpz(1))
  zn = S(n*dn)
  zd = S(d*dd)
  cn = content(zn)
  cd = content(zd)
  zn = divexact(zn, cn)
  zd = divexact(zd, cd)
  cn *= dd
  cd *= dn
  g = gcd(cn, cd)
  cn = divexact(cn, g)
  cd = divexact(cd, g)
  return cn*zn, cd*zd
end

function (S::FmpzPolyRing)(a::Generic.Rat{fmpq})
  n, d = integral_split(a, S)
  @assert isone(d)
  return n
end

end

using .HessMain

#=
  this should work:

Hecke.example("Round2.jl")

?GenericRound2

Qt, t = RationalFunctionField(QQ, "t")
Qtx, x = PolynomialRing(Qt, "x")
F, a = FunctionField(x^6+27*t^2+108*t+108, "a")
integral_closure(parent(denominator(t)), F)
integral_closure(Localization(Qt, degree), F)
integral_closure(Hecke.Globals.Zx, F)
basis(ans, F)
derivative(F.pol)(gen(F)) .* ans #should be integral

k, a = wildanger_field(3, 8*13)
integral_closure(ZZ, k)
integral_closure(Localization(ZZ, 2), k)

more interesting and MUCH harder:

G, b = FunctionField(x^6 + (140*t - 70)*x^3 + 8788*t^2 - 8788*t + 2197, "b")

=#


text:
# "Send and receive" is an integration test that verifies that the client can send and receive
# messages.
#
# Five messages are sent from the client to the server, and five messages are sent from the server
# to the client.

include("Stubs.jl")

using DandelionWebSockets
using DandelionWebSockets: HTTPHandshake
using .Stubs
using .Stubs: write, read, InProcessIO
using Test
using Random

# These are the actions that the server will take.
serverscript = [
    WaitForOpen(),
    SendTextFrame("1234"),
    ShortWait(),
    SendTextFrame("2345"),
    ShortWait(),
    SendTextFrame("5678"),
]

clientscript = [
    WaitForOpen(),
    SendTextFrame("abcd"),
    ShortWait(),
    SendTextFrame("bcde"),
    ShortWait(),
    SendTextFrame("cdef"),
    ShortWait(),
    SendTextFrame("defg"),
    ShortWait(),
    SendTextFrame("efgh"),
    ShortWait(),
    CloseConnection()
]

# Create server, and the handshake adapter that creates the connection when the client requests a
# handshake.
server = ScriptedServer(serverscript)
handshakeadapter = InProcessHandshakeAdapter(server)

# Create the WebSocket client, with a custom handshake adapter.
# Instead of making an actual HTTP request, the custom handshake adapter creates an in-process
# `IO` that is used in place of an actual network socket.
# Note that the HTTP handshake logic is still the same, it's only the actual transport method that
# is different.
client = WSClient(; handshake=HTTPHandshake(RandomDevice(), handshakeadapter))

# Create a client that runs the script above.
clienthandler = ScriptedClientHandler(clientscript)

# Request that the connection be opened, which starts the scripts for both server and client.
wsconnect(client, "ws://the/uri/does/not/matter/here", clienthandler)

# Wait for both client and server to close.
waitforscriptdone(server, clienthandler)

@testset "Sent messages    " begin
    @test server.statistics.texts_received == 5
end

@testset "Received messages" begin
    @test clienthandler.statistics.texts_received == 3
end


text:
module NonconvexUtils

export  ForwardDiffFunction,
        AbstractDiffFunction,
        AD,
        TraceFunction,
        CustomGradFunction,
        LazyJacobian,
        CustomHessianFunction,
        ImplicitFunction

using ChainRulesCore, AbstractDifferentiation, ForwardDiff, LinearAlgebra
using Zygote, LinearMaps, IterativeSolvers, NonconvexCore
using NonconvexCore: flatten
using MacroTools

struct AbstractDiffFunction{F, B} <: Function
    f::F
    backend::B
end
ForwardDiffFunction(f) = AbstractDiffFunction(f, AD.ForwardDiffBackend())
(f::AbstractDiffFunction)(x) = f.f(x)
function ChainRulesCore.rrule(
    f::AbstractDiffFunction, x::AbstractVector,
)
    v, (∇,) = AbstractDifferentiation.value_and_gradient(f.backend, f.f, x)
    return v, Δ -> (NoTangent(), Δ * ∇)
end

struct TraceFunction{F, V} <: Function
    f::F
    trace::V
    on_call::Bool
    on_grad::Bool
end
function TraceFunction(f; on_call::Union{Bool, Nothing} = nothing, on_grad::Union{Bool, Nothing} = nothing)
    if on_call === on_grad === nothing
        _on_call = true
        _on_grad = true
    elseif on_call === nothing
        _on_call = !on_grad
        _on_grad = on_grad
    elseif on_grad === nothing
        _on_call = on_call
        _on_grad = !on_call
    else
        _on_call = on_call
        _on_grad = on_grad
    end
    return TraceFunction(f, Any[], _on_call, _on_grad)
end
function (tf::TraceFunction)(x)
    v = tf.f(x)
    if tf.on_call
        push!(tf.trace, (input = copy(x), output = copy(v)))
    end
    return v
end
function ChainRulesCore.rrule(rc::RuleConfig, tf::TraceFunction, x)
    v, pb = ChainRulesCore.rrule_via_ad(rc, tf.f, x)
    return v, Δ -> begin
        Δin = pb(Δ)
        g = Δin[2] isa Array ? Δin[2] : Δin[2].val.f()
        if tf.on_grad
            push!(tf.trace, (input = copy(x), output = copy(v), grad = copy(g)))
        end
        return (Δin[1], g)
    end
end

struct LazyJacobian{symmetric, J1, J2}
    jvp::J1
    jtvp::J2
end
function LazyJacobian(; jvp=nothing, jtvp=nothing, symmetric=false)
    return LazyJacobian{symmetric}(jvp, jtvp)
end
function LazyJacobian{symmetric}(jvp = nothing, jtvp = nothing) where {symmetric}
    if jvp === jtvp === nothing
        throw(ArgumentError("Both the jvp and jtvp operators cannot be nothing."))
    end
    if symmetric
        if jvp !== nothing
            _jtvp = _jvp = jvp
        else
            _jvp = _jtvp = jtvp
        end
    else
        _jvp = jvp
        _jtvp = jtvp
    end
    return LazyJacobian{symmetric, typeof(_jvp), typeof(_jtvp)}(_jvp, _jtvp)
end

struct LazyJacobianTransposed{J}
    j::J
end

LinearAlgebra.adjoint(j::LazyJacobian{false}) = LazyJacobianTransposed(j)
LinearAlgebra.transpose(j::LazyJacobian{false}) = LazyJacobianTransposed(j)
LinearAlgebra.adjoint(j::LazyJacobian{true}) = j
LinearAlgebra.transpose(j::LazyJacobian{true}) = j
LinearAlgebra.adjoint(j::LazyJacobianTransposed) = j.j
LinearAlgebra.transpose(j::LazyJacobianTransposed) = j.j

LinearAlgebra.:*(j::LazyJacobian, v::AbstractVecOrMat) = j.jvp(v)
LinearAlgebra.:*(v::AbstractVecOrMat, j::LazyJacobian) = j.jtvp(v')'
LinearAlgebra.:*(j::LazyJacobianTransposed, v::AbstractVecOrMat) = (v' * j')'
LinearAlgebra.:*(v::AbstractVecOrMat, j::LazyJacobianTransposed) = (j' * v')'

struct CustomGradFunction{F, G} <: Function
    f::F
    g::G
end
(f::CustomGradFunction)(x) = f.f(x)
function ChainRulesCore.rrule(f::CustomGradFunction, x)
    return f.f(x), Δ -> begin
        G = f.g(x)
        if G isa AbstractVector
            return (NoTangent(), G * Δ)
        else
            return (NoTangent(), G' * Δ)
        end
    end
end

struct CustomHessianFunction{F, G, H} <: Function
    f::F
    g::G
    h::H
    function CustomHessianFunction(
        f::F, g::G, h::H; hvp = false,
    ) where {F, G, H}
        _h = hvp ? x -> LazyJacobian{true}(v -> h(x, v)) : h
        return new{F, G, typeof(_h)}(f, g, _h)
    end
end
(to::CustomHessianFunction)(x) = to.f(x)
function ChainRulesCore.rrule(f::CustomHessianFunction, x)
    g = CustomGradFunction(f.g, f.h)
    return f.f(x), Δ -> (NoTangent(), g(x) * Δ)
end

# Parameters x, variables y, residuals f
struct ImplicitFunction{matrixfree, F, C, L, T} <: Function
    # A function which takes x as input and returns a tuple (ystar, df/dy) such that f(x, ystar) = 0. df/dy is optional and can be replaced by nothing to compute it via automatic differentiation. Jacobian should only be returned if it's more cheaply available than using AD, e.g. when using BFGS approximation of the Hessian in IPOPT.
    forward::F
    # The conditions function f(x, y) which must be 0 at ystar. Note that variables which don't show up in x and are closed over instead will be assumed to have no effect on the optimal solution. So it's the user's responsibility to ensure x includes all the interesting variables to be differentiated with respect to.
    conditions::C
    # A linear system solver to solve df/dy' \ v
    linear_solver::L
    # The acceptable tolerance for f(x, ystar) to use the implicit function theorem at x
    tol::T
    # A booolean to decide whether or not to error if the tolerance is violated, i.e. norm(f(x, ystar)) > tol. If false, we return a gradient of NaNs.
    error_on_tol_violation::Bool
end
function ImplicitFunction(
    forward::F, conditions::C; tol::T = 1e-5, error_on_tol_violation = false, matrixfree = false, linear_solver::L = _default_solver(matrixfree),
) where {F, C, L, T}
    return ImplicitFunction{matrixfree, F, C, L, T}(
        forward, conditions, linear_solver, tol, error_on_tol_violation,
    )
end

function _default_solver(matrixfree)
    if matrixfree
        return (A, b) -> begin
            L = LinearMap(A, length(b))
            return gmres(L, b)
        end
    else
        return (A, b) -> A \ b
    end
end

(f::ImplicitFunction)(x) = f.forward(x)[1]
(f::ImplicitFunction)() = f.forward()[1]

function ChainRulesCore.rrule(
    rc::RuleConfig, f::ImplicitFunction{matrixfree}, x,
) where {matrixfree}
    ystar, _dfdy = f.forward(x)
    flat_ystar, unflatten_y = flatten(ystar)
    forward_returns_jacobian = _dfdy !== nothing
    if forward_returns_jacobian
        dfdy = _dfdy
        if matrixfree
            # y assumed flat if dfdy is passed in
            pby = v -> dfdy' * v
        else
            pby = nothing
        end
    else
        _conditions_y = flat_y -> begin
            return flatten(f.conditions(x, unflatten_y(flat_y)))[1]
        end
        if matrixfree
            dfdy = nothing
            _, _pby = rrule_via_ad(rc, _conditions_y, flat_ystar)
            pby = v -> _pby(v)[2]
        else
            # Change this to AbstractDifferentiation
            dfdy = Zygote.jacobian(_conditions_y, flat_ystar)[1]
            pby = nothing
        end
    end
    _conditions_x = (conditions, x) -> begin
        return flatten(conditions(x, ystar))[1]
    end
    residual, pbx = rrule_via_ad(rc, _conditions_x, f.conditions, x)
    return ystar, ∇ -> begin
        if norm(residual) > f.tol && f.error_on_tol_violation
            throw(ArgumentError("The acceptable tolerance for the implicit function theorem is not satisfied for the current problem. Please double check your function definition, increase the tolerance, or set `error_on_tol_violation` to false to ignore the violation and return `NaN`s for the gradient."))
        end
        if matrixfree
            ∇f, ∇x = Base.tail(pbx(f.linear_solver(pby, -flatten(∇)[1])))
        else
            ∇f, ∇x = Base.tail(pbx(f.linear_solver(dfdy', -flatten(∇)[1])))
        end
        ∇imf = Tangent{typeof(f)}(
            conditions = Tangent{typeof(f.conditions)}(;
                ChainRulesCore.backing(∇f)...,
            ),
        )
        if norm(residual) <= f.tol
            return (∇imf, ∇x)
        else
            return (nanlike(∇imf), nanlike(∇x))
        end
    end
end

function ChainRulesCore.rrule(
    rc::RuleConfig, f::ImplicitFunction{matrixfree},
) where {matrixfree}
    ystar, _dfdy = f.forward()
    flat_ystar, unflatten_y = flatten(ystar)
    forward_returns_jacobian = _dfdy !== nothing
    if forward_returns_jacobian
        dfdy = _dfdy
        if matrixfree
            # y assumed flat if dfdy is passed in
            pby = v -> dfdy' * v
        else
            pby = nothing
        end
    else
        _conditions_y = flat_y -> begin
            return flatten(f.conditions(unflatten_y(flat_y)))[1]
        end
        if matrixfree
            dfdy = nothing
            _, _pby = rrule_via_ad(rc, _conditions_y, flat_ystar)
            pby = v -> _pby(v)[2]
        else
            # Change this to AbstractDifferentiation
            dfdy = Zygote.jacobian(_conditions_y, flat_ystar)[1]
            pby = nothing
        end
    end
    _conditions = (conditions) -> begin
        return flatten(conditions(ystar))[1]
    end
    residual, pbf = rrule_via_ad(rc, _conditions, f.conditions)
    return ystar, ∇ -> begin
        if norm(residual) > f.tol && f.error_on_tol_violation
            throw(ArgumentError("The acceptable tolerance for the implicit function theorem is not satisfied for the current problem. Please double check your function definition, increase the tolerance, or set `error_on_tol_violation` to false to ignore the violation and return `NaN`s for the gradient."))
        end
        if matrixfree
            ∇f = pbf(f.linear_solver(pby, -flatten(∇)[1]))[2]
        else
            ∇f = pbf(f.linear_solver(dfdy', -flatten(∇)[1]))[2]
        end
        ∇imf = Tangent{typeof(f)}(
            conditions = Tangent{typeof(f.conditions)}(;
                ChainRulesCore.backing(∇f)...,
            ),
        )
        if norm(residual) <= f.tol
            return (∇imf,)
        else
            return (nanlike(∇imf),)
        end
    end
end

function nanlike(x)
    flat, un = flatten(x)
    return un(similar(flat) .= NaN)
end

macro ForwardDiff_frule(sig)
    _fd_frule(sig)
end
function _fd_frule(sig)
    MacroTools.@capture(sig, f_(x__))
    return quote
        function $(esc(f))($(esc.(x)...))
            f = $(esc(f))
            x = ($(esc.(x)...),)
            flatx, unflattenx = flatten(x)
            CS = length(ForwardDiff.partials(first(flatx)))
            flat_xprimals = ForwardDiff.value.(flatx)
            flat_xpartials = reduce(vcat, transpose.(ForwardDiff.partials.(flatx)))

            xprimals = unflattenx(flat_xprimals)
            xpartials1 = unflattenx(flat_xpartials[:,1])

            yprimals, ypartials1 = ChainRulesCore.frule(
                (NoTangent(), xpartials1...), f, xprimals...,
            )
            flat_yprimals, unflatteny = flatten(yprimals)
            flat_ypartials1, _ = flatten(ypartials1)
            flat_ypartials = hcat(reshape(flat_ypartials1, :, 1), ntuple(Val(CS - 1)) do i
                xpartialsi = unflattenx(flat_xpartials[:, i+1])
                _, ypartialsi = ChainRulesCore.frule((NoTangent(), xpartialsi...), f, xprimals...)
                return flatten(ypartialsi)[1]
            end...)

            T = ForwardDiff.tagtype(eltype(flatx))
            flaty = ForwardDiff.Dual{T}.(
                flat_yprimals, ForwardDiff.Partials.(NTuple{CS}.(eachrow(flat_ypartials))),
            )
            return unflatteny(flaty)
        end
    end
end

end


text:
# by Lazaro Alonso
using GLMakie, GeoMakie, CairoMakie, Downloads
using GeoJSON, GeoInterface
CairoMakie.activate!()
let
    # data from
    # https://github.com/telegeography/www.submarinecablemap.com
    urlPoints = "https://raw.githubusercontent.com/telegeography/www.submarinecablemap.com/master/web/public/api/v3/landing-point/landing-point-geo.json"
    urlCables = "https://raw.githubusercontent.com/telegeography/www.submarinecablemap.com/master/web/public/api/v3/cable/cable-geo.json"

    landPoints = Downloads.download(urlPoints, IOBuffer())
    landCables = Downloads.download(urlCables, IOBuffer())

    land_geoPoints = GeoJSON.read(seekstart(landPoints))
    land_geoCables = GeoJSON.read(seekstart(landCables))

    toPoints = GeoMakie.geo2basic(land_geoPoints)
    #toLines = GeoMakie.geo2basic(land_geoCables) # this should probably be supported.
    feat = GeoInterface.features(land_geoCables)
    toLines = GeoInterface.coordinates.(GeoInterface.geometry.(feat))
    splitLines = []
    for i in 1:length(toLines)
        for j in 1:length(toLines[i])
            push!(splitLines, Point2f0.(toLines[i][j]))
        end
    end
    # broken lines at -180 and 180... they should
    # be the same line and be in the same array.

    function plotCables()
        fig, ax, = scatter(toPoints, markersize= 5, color = 1:length(toPoints),
            colormap = :plasma, figure = (;resolution = (1200,700), fontsize = 24))
        [lines!(ax,splitLines[i], linewidth=0.85) for i in 1:length(splitLines)]
        lines!(GeoMakie.coastlines(), color = :white, linewidth = 0.35)
        limits!(ax,-185,185,-95,95)
        ax.xticks = -180:60:180
        ax.yticks = -90:30:90
        ax.xtickformat = "{:d}ᵒ"
        ax.ytickformat = "{:d}ᵒ"
        fig
    end
    fig = with_theme(plotCables, theme_dark())
    save(joinpath(@__DIR__, "output", "submarineCables.png"), fig) # HIDE
end
using Pkg # HIDE
Pkg.status(["GLMakie", "GeoMakie", "CairoMakie", "Downloads", "GeoJSON", "GeoInterface"]) # HIDE

max_stars_repo_name: MisterBiggs/BeautifulMakie

max_issues_repo_name: MisterBiggs/BeautifulMakie

max_forks_repo_name: MisterBiggs/BeautifulMakie


text:
using FileIO
using FactCheck
using Compat

if VERSION < v"0.4.0-dev"
    import FileIO.Pair
end

context("OS") do
    if is_linux()
        @fact FileIO.applies_to_os(FileIO.Linux)   --> true
        @fact FileIO.applies_to_os(FileIO.OSX)     --> false
        @fact FileIO.applies_to_os(FileIO.Unix)    --> true
        @fact FileIO.applies_to_os(FileIO.Windows) --> false
    end
    if is_apple()
        @fact FileIO.applies_to_os(FileIO.Linux)   --> false
        @fact FileIO.applies_to_os(FileIO.OSX)     --> true
        @fact FileIO.applies_to_os(FileIO.Unix)    --> true
        @fact FileIO.applies_to_os(FileIO.Windows) --> false
    end
    if is_windows()
        @fact FileIO.applies_to_os(FileIO.Linux)   --> false
        @fact FileIO.applies_to_os(FileIO.OSX)     --> false
        @fact FileIO.applies_to_os(FileIO.Unix)    --> false
        @fact FileIO.applies_to_os(FileIO.Windows) --> true
    end
end

# Before we bork things, make a copy
ext2sym = copy(FileIO.ext2sym)
magic_list = copy(FileIO.magic_list)
sym2info = copy(FileIO.sym2info)

module LoadTest1
import FileIO: @format_str, File
load(file::File{format"MultiLib"}) = error()

save(file::File{format"MultiLib"}) = open(file, "w") do s
    write(s, magic(format"MultiLib"))  # Write the magic bytes
    write(s, 0)
end

end
module LoadTest2
import FileIO: @format_str, File, magic
load(file::File{format"MultiLib"}) = 42

save(file::File{format"MultiLib"}) = open(file, "w") do s
    write(s, magic(format"MultiLib"))  # Write the magic bytes
    write(s, 42)
end

end

try
    empty!(FileIO.ext2sym)
    empty!(FileIO.magic_list)
    empty!(FileIO.sym2info)

    context("DataFormat") do
        @fact DataFormat{:CSV} --> format"CSV"
        @fact unknown(format"CSV") --> false
        @fact unknown(format"UNKNOWN") --> true

        add_format(format"CSV", UInt8[], ".csv")
        @fact info(format"CSV") --> ((), ".csv")
        add_format(format"FOO", (), ".foo")  # issue #17
        @fact_throws info(format"OOPS")
        @fact FileIO.ext2sym[".csv"] --> :CSV
        del_format(format"FOO")
        @fact FileIO.magic_list --> [Pair((), :CSV)]
        del_format(format"CSV")
        @fact isempty(FileIO.ext2sym) --> true
        @fact isempty(FileIO.magic_list) --> true
        @fact isempty(FileIO.sym2info) --> true
        @fact_throws info(format"CSV")

        add_format(format"JUNK", "JUNK", [".jnk",".junk",".JNK"])
        @fact info(format"JUNK") --> (tuple(b"JUNK"...), [".jnk",".junk",".JNK"])
        @fact FileIO.ext2sym[".jnk"]  --> :JUNK
        @fact FileIO.ext2sym[".junk"] --> :JUNK
        @fact FileIO.ext2sym[".JNK"]  --> :JUNK
        @fact FileIO.magic_list --> [Pair((0x4a,0x55,0x4e,0x4b),:JUNK)]
        @fact_throws add_format(format"JUNK2", "JUNK", ".jnk2")  # magic bytes already registered

    end

    context("streams") do
        io = IOBuffer()
        s = Stream(format"JUNK", io)
        @fact typeof(s) --> Stream{DataFormat{:JUNK}, IOBuffer}
        @fact isnull(filename(s)) --> true
        @fact_throws FileIO.file!(s)
        s = Stream(format"JUNK", io, "junk.jnk")
        @fact get(filename(s)) --> "junk.jnk"
        s = Stream(format"JUNK", io, Nullable("junk2.jnk"))
        @fact get(filename(s)) --> "junk2.jnk"
    end

    context("query") do
        # Streams
        io = IOBuffer()
        write(io, "Weird format")
        seek(io, 0)
        q = query(io)
        @fact unknown(q) --> true

        # Short "file"
        truncate(io, 0)
        write(io, "S")
        seek(io, 0)
        q = query(io)
        @fact unknown(q) --> true

        truncate(io, 0)
        write(io, "JUNK and some more stuff")
        seek(io, 0)
        q = query(io)
        @fact typeof(q) --> Stream{format"JUNK",typeof(io)}
        @fact unknown(q) --> false
        @fact isnull(file_extension(q)) --> true

        # File with correct extension
        str = takebuf_string(io)
        fn = string(tempname(), ".jnk")
        open(fn, "w") do file
            write(file, str)
        end
        q = query(fn)
        @fact typeof(q) --> File{format"JUNK"}
        @fact file_extension(q) --> ".jnk"

        rm(fn)

        # File with erroneous extension
        fn = string(tempname(), ".csv")
        open(fn, "w") do file
            write(file, str)
        end
        q = query(fn)
        @fact typeof(q) --> File{format"JUNK"}
        @fact file_extension(q) --> ".csv"
        rm(fn)

        # Format with no magic bytes
        add_format(format"BAD", (), ".bad")
        fn = string(tempname(), ".bad")
        open(fn, "w") do file
            write(file, "Here's some data")
        end
        q = query(fn)
        @fact typeof(q) --> File{format"BAD"}
        @fact file_extension(q) --> ".bad"
        rm(fn)

        q = query( "some_non_existant_file.bad")
        @fact typeof(q) --> File{format"BAD"}

        # Unknown extension
        fn = string("tempname", ".wrd")
        open(fn, "w") do file
            write(file, "More data")
        end
        @fact unknown(query(fn)) --> true
        rm(fn)

        add_format(format"DOUBLE_1", "test1", ".double")
        add_format(format"DOUBLE_2", "test2", ".double")

        fn = string(tempname(), ".double")
        open(fn, "w") do file
            write(file, "test1")
        end
        q = query(fn)
        @fact typeof(q) --> File{format"DOUBLE_1"}
        rm(fn)

        add_format(format"MAGIC", "this so magic", ".mmm")
        q = query( "some_non_existant_file.mmm")
        @fact typeof(q) --> File{format"MAGIC"}

        add_format(format"DOUBLE_MAGIC", (UInt8[0x4d,0x4d,0x00,0x2a], UInt8[0x4d,0x4d,0x00]), ".dd2")

        fn = string(tempname(), ".dd2")
        open(fn, "w") do file
            write(file, UInt8[0x4d,0x4d,0x00,0x2a])
            write(file, randstring(19))
        end
        q = query(fn)
        @fact typeof(q) --> File{format"DOUBLE_MAGIC"}
        io = open(q)
        skipmagic(io)
        @fact position(io) --> 4
        close(io)
        rm(fn)

        open(fn, "w") do file
            write(file, UInt8[0x4d,0x4d,0x00])
            write(file, randstring(19))
        end
        q = query(fn)
        @fact typeof(q) --> File{format"DOUBLE_MAGIC"}
        io = open(q)
        @fact file_extension(q) --> ".dd2"
        skipmagic(io)
        @fact position(io) --> 3
        close(io)
        open(fn, "w") do file
            write(file, randstring(19)) # corrupt magic bytes
        end
        open(fn, "r") do file
            @fact_throws skipmagic(file)
        end
        rm(fn)
        lene0 = length(FileIO.ext2sym)
        lenm0 = length(FileIO.magic_list)
        del_format(format"DOUBLE_MAGIC")
        @fact lene0 - 1 --> length(FileIO.ext2sym)
        @fact lenm0 - 2 --> length(FileIO.magic_list)
    end

    del_format(format"JUNK")  # This triggers del_extension for multiple extensions

    context("multiple libs") do
        lensave0 = length(FileIO.sym2saver)
        lenload0 = length(FileIO.sym2loader)
        OSKey = is_apple() ? FileIO.OSX : is_windows() ? FileIO.Windows : is_linux() ? FileIO.Linux : error("os not supported")
        add_format(
            format"MultiLib",
            UInt8[0x42,0x4d],
            ".mlb",
            [:LoadTest1, FileIO.LOAD, OSKey],
            [:LoadTest2]
        )
        stderr_copy = STDERR
        rd,wr = redirect_stderr()
        @fact lensave0 + 1 --> length(FileIO.sym2saver)
        @fact lenload0 + 1 --> length(FileIO.sym2loader)
        @fact length(FileIO.sym2loader[:MultiLib]) --> 2
        @fact length(FileIO.sym2saver[:MultiLib]) --> 1
        fn = string(tempname(), ".mlb")
        save(fn)
        x = load(fn)
        open(query(fn), "r") do io
            skipmagic(io)
            a = read(io, Int)
            @fact a --> 42 #make sure that LoadTest2 is used for saving, even though its at position 2
        end
        @fact isdefined(:LoadTest1) --> true # first module should load first but fail
        @fact x --> 42
        rm(fn)
        close(rd);close(wr)
        redirect_stderr(stderr_copy)
    end
finally
    # Restore the registry
    empty!(FileIO.ext2sym)
    empty!(FileIO.magic_list)
    empty!(FileIO.sym2info)

    merge!(FileIO.ext2sym, ext2sym)
    append!(FileIO.magic_list, magic_list)
    merge!(FileIO.sym2info, sym2info)
end

file_dir = joinpath(dirname(@__FILE__), "files")
context("STL detection") do
    q = query(joinpath(file_dir, "ascii.stl"))
    @fact typeof(q) --> File{format"STL_ASCII"}
    q = query(joinpath(file_dir, "binary_stl_from_solidworks.STL"))
    @fact typeof(q) --> File{format"STL_BINARY"}
    open(q) do io
        @fact position(io) --> 0
        skipmagic(io)
        @fact position(io) --> 0 # no skipping for functions
    end
end
context("PLY detection") do
    q = query(joinpath(file_dir, "ascii.ply"))
    @fact typeof(q) --> File{format"PLY_ASCII"}
    q = query(joinpath(file_dir, "binary.ply"))
    @fact typeof(q) --> File{format"PLY_BINARY"}

end
context("Multiple Magic bytes") do
    q = query(joinpath(file_dir, "magic1.tiff"))
    @fact typeof(q) --> File{format"TIFF"}
    q = query(joinpath(file_dir, "magic2.tiff"))
    @fact typeof(q) --> File{format"TIFF"}
    open(q) do io
        @fact position(io) --> 0
        skipmagic(io)
        @fact position(io) --> 4
    end
end
context("AVI Detection") do
    open(joinpath(file_dir, "bees.avi")) do s
        @fact FileIO.detectavi(s) --> true
    end
    open(joinpath(file_dir, "sin.wav")) do s
        @fact FileIO.detectavi(s) --> false
    end
    open(joinpath(file_dir, "magic1.tiff")) do s
        @fact FileIO.detectavi(s) --> false
    end
    q = query(joinpath(file_dir, "bees.avi"))
    @fact typeof(q) --> File{format"AVI"}
end
context("RDA detection") do
    q = query(joinpath(file_dir, "minimal_ascii.rda"))
    @fact typeof(q) --> File{format"RData"}
    open(q) do io
        @fact position(io) --> 0
        @fact FileIO.detect_rdata(io) --> true
        @fact position(io) --> 5
    end
end
context("Format with function for magic bytes") do
    add_format(format"FUNCTION_FOR_MAGIC_BYTES", x -> 0x00, ".wav", [:WAV])
    del_format(format"FUNCTION_FOR_MAGIC_BYTES")
end


text:
module DashReactFontAwesomeIconPicker
using Dash

const resources_path = realpath(joinpath( @__DIR__, "..", "deps"))
const version = "0.0.1"

include("jl/dashiconpicker.jl")
include("jl/dashiconpickeritem.jl")

function __init__()
    DashBase.register_package(
        DashBase.ResourcePkg(
            "dash_react_font_awesome_icon_picker",
            resources_path,
            version = version,
            [
                DashBase.Resource(
    relative_package_path = "dash_react_font_awesome_icon_picker.min.js",
    external_url = nothing,
    dynamic = nothing,
    async = nothing,
    type = :js
),
DashBase.Resource(
    relative_package_path = "dash_react_font_awesome_icon_picker.min.js.map",
    external_url = nothing,
    dynamic = true,
    async = nothing,
    type = :js
)
            ]
        )

    )
end
end