text:
# ==============================================================================
using PyPlot, AcousticWaveLab
nz = 201; nx = 373; ext = 30; iflag = 2; dz = 10.; dx = 10.; dt = 1e-3;
v  = 3000. * ones(nz, nx);
v[101:end, :] = 4000.
# initialize source term
pos = [1 180]; f0 = 30.0; ot = [0.0];
src = InitSources(pos, f0, ot, dt);

# discretize spatial derivative operator
vmax = maximum(v); vmin = minimum(v);
fidMtx = InitFidMtx(nz, nx, ext, iflag, dz, dx, dt, vmax, vmin, f0, v);

ix = collect(1:5:nx); iz = 1 * ones(Int64, length(ix));
pos = hcat(iz, ix)
shot = MultiStepForward(pos, src[1], fidMtx, print_flag=true)
SeisPlot(shot.d, pclip=98)
# ==============================================================================
# using PyPlot, AcousticWaveLab
# nz = 201; nx = 201; ext = 20; iflag = 1; dz = 10.; dx = 10.; dt = 2e-3;
# v  = 3000. * ones(nz, nx);
# # initialize source term
# pos = [101 101]; f0 = 30.0; ot = [0.0];
# src = InitSources(pos, f0, ot, dt);
# # discretize spatial derivative operator
# vmax = maximum(v); vmin = minimum(v);
# fidMtx = InitFidMtx(nz, nx, ext, iflag, dz, dx, dt, vmax, vmin, f0, v);
# path = "/Users/wenlei/Desktop/acoustic.bin"
# MultiStepForward(path, src, fidMtx, tmax=1.0);
#
# iz = 5 * ones(Int64, 201); ix = collect(1:201);
# pos = hcat(iz, ix)
# shot = fullWfd2Shot(pos, path);
# shot1 = MultiStepForward(pos, src[1], fidMtx, tmax=1.0)
#
# path = "/Users/wenlei/Desktop/adjoint.bin"
# MultiStepAdjoint!(path, shot, fidMtx)
#
#
# using PyPlot
# imshow(read(path, 500))
#
# iz = 5 * ones(Int64, 201); ix = collect(1:201);
# pos = hcat(iz, ix)
# shot = wfd2shot(pos, path);
# figure()
# plot(shot.d[:,71])

# ==============================================================================
# dot product test for one step
# spt1 = InitSnapShot(101, 101, nz, nx, ext, 1, dt, 1);
# spt1.vz = randn(length(spt1.vz));
# spt1.vx = randn(length(spt1.vx));
# spt1.pz = randn(length(spt1.pz));
# spt1.px = randn(length(spt1.px));
# spt2 = InitSnapShot(101, 101, nz, nx, ext, 1, dt, 1);
# OneStepForward!(spt2, spt1, fidMtx);
#
# tmp2 = InitSnapShot(101, 101, nz, nx, ext, 1, dt, 2)
# tmp2.vz = randn(length(tmp2.vz))
# tmp2.vx = randn(length(tmp2.vx))
# tmp2.pz = randn(length(tmp2.pz))
# tmp2.px = randn(length(tmp2.px))
# tmp1 = InitSnapShot(101, 101, nz, nx, ext, 1, dt, 2)
# OneStepAdjoint!(tmp1, tmp2, fidMtx)
#
# x  = vcat(spt1.vz, spt1.vx, spt1.pz, spt1.px)
# x1 = vcat(tmp1.vz, tmp1.vx, tmp1.pz, tmp1.px)
# y  = vcat(spt2.vz, spt2.vx, spt2.pz, spt2.px)
# y1 = vcat(tmp2.vz, tmp2.vx, tmp2.pz, tmp2.px)
#
# x'*x1 - y'*y1

# ==============================================================================
# nz = 201; nx = 201; dt = 2e-3; iflag = 1; ext = 20
# if iflag ==1
#    zupper = ext
#    Nz = nz + 2*ext
# elseif iflag == 2
#    Nz = nx +   ext
#    zupper = ext
# end
# Nx = nx + 2*ext
# spt1 = InitSnapShot(101, 101, nz, nx, ext, iflag, dt, 1);
# spt1.vz = randn(length(spt1.vz));
# spt1.vx = randn(length(spt1.vx));
# spt1.pz = randn(length(spt1.pz));
# spt1.px = randn(length(spt1.px));
#
# iz = collect(1:5:201); ix = collect(1:3:201);
# nrz= length(iz); nrx = length(ix);
# iz = vec(repmat(iz, nrx, 1));
# ix = vec(repmat(ix', (nrz), 1));
# pos = hcat(iz, ix);
# shot1 = InitShot(101, 101, pos, 0.0, 1, dt)
# spt2shot!(shot1, spt1, Nz, Nx, ext, zupper)
#
#
# shot2 = InitShot(101, 101, pos, 0.0, 1, dt)
# shot2.d = randn(1, size(shot2.d, 2))
# spt2 = InitSnapShot(101, 101, nz, nx, ext, iflag, dt, 1);
# AddShot2SnapShot!(spt2, shot2, Nz, Nx, ext, zupper)
#
# x  = vcat(spt1.vz, spt1.vx, spt1.pz, spt1.px)
# x1 = vcat(spt2.vz, spt2.vx, spt2.pz, spt2.px)
# y  = vec(shot1.d)
# y1 = vec(shot2.d)
# x'*x1 - y'*y1

# ==============================================================================
# using AcousticWave
# nz = 201; nx = 201; ext = 20; iflag = 1; dz = 10.; dx = 10.; dt = 2e-3; f0 = 30.0
# v  = 3000. * ones(nz, nx);
# vmax = maximum(v); vmin = minimum(v);
# fidMtx = InitFidMtx(nz, nx, ext, iflag, dz, dx, dt, vmax, vmin, f0, v);
# # initialize source term
# ix  = collect(1:2:201); iz = 101*ones(Int64, length(ix));
# pos = hcat(iz, ix)
# f0 = 30.0; ot = zeros(size(pos,1));
# src = InitSources(pos, f0, ot, dt);
# path = "/Users/wenlei/Desktop/src2spts.bin"
# src2SnapShots(path, src, nz, nx, ext, iflag);
#
# # position of receiver
# irz = []; irx = []
# for ix = 1 : nx
#     for iz = 1 : nz
#         if rand() <= 0.1
#            irz = push!(irz, iz)
#            irx = push!(irx, ix)
#         end
#     end
# end
# pos_rec = convert(Array{Int64,2}, hcat(irz, irx))
# shot  = MultiStepForward(pos_rec, path, fidMtx, tmax=1.0)
# path_wfd = "/Users/wenlei/Desktop/acoustic.bin"
# MultiStepForward(path_wfd, src , fidMtx, tmax=1.0)
# shot1 = fullWfd2Shot(pos_rec, path_wfd)

# ==============================================================================
# using AcousticWave
# nz = 201; nx = 201; ext = 20; iflag = 1; dz = 10.; dx = 10.; dt = 2e-3; f0 = 30.0
# v  = 3000. * ones(nz, nx);
# vmax = maximum(v); vmin = minimum(v);
# fidMtx = InitFidMtx(nz, nx, ext, iflag, dz, dx, dt, vmax, vmin, f0, v);
# # initialize source term
# isz = []; isx = []
# for ix = 1 : nx
#     for iz = 1 : nz
#         if rand() <= 0.1
#            isz = push!(isz, iz)
#            isx = push!(isx, ix)
#         end
#     end
# end
# pos_src = convert(Array{Int64,2}, hcat(isz, isx))
# f0 = 30.0; ot = zeros(size(pos_src,1));
# src = InitSources(pos_src, f0, ot, dt);
# tmax = 1; nt = round(Int64, tmax/dt+1);
# w = Ricker(f0, dt); ns = nt + 1 - length(w)
# for isrc = 1 : length(src)
#     src[isrc].nt= nt
#     src[isrc].p = conv(randn(ns),w)
# end
# path = "/Users/wenlei/Desktop/src2spts.bin"
# src2SnapShots(path, src, nz, nx, ext, iflag);
# # position of receiver
# irx = collect(1:2:nx); irz = 5*ones(Int64,length(irx));
# pos_rec = hcat(irz, irx)
# shot  = MultiStepForward(pos_rec, path, fidMtx, tmax=1.0)
#
#
# path_adj = "/Users/wenlei/Desktop/adjoint.bin"
# shot1 = InitShot(0, 0, pos_rec, 0.0, nt, dt)
# for irec = 1 : length(shot1.irz)
#     shot1.d[:,irec] = conv(randn(ns), w)
# end
# MultiStepAdjoint_test(path_adj, shot1, fidMtx)
#
# L = 0.0
# for ir = 1 : length(shot.irx)
#     L = L + dot(shot.d[:,ir], shot1.d[:,ir])
# end
#
# T = 0.0
# for it = 1 : nt
#     spt1 = readSnapShot(path    , it)
#     spt2 = readSnapShot(path_adj, it)
#     T = T + dot(spt1.vz, spt2.vz) + dot(spt1.vx, spt2.vx) + dot(spt1.pz, spt2.pz) + dot(spt1.px, spt2.px)
# end

# ==============================================================================
# using AcousticWave
# path = "/Users/wenlei/Desktop/acoustic.bin"
# nz = 201; nx = 201; ext = 20; iflag = 1; dz = 10.; dx = 10.; dt = 2e-3;
# v  = 3000. * ones(nz, nx);
# # initialize source term
# pos = [101 101]; f0 = 30.0; ot = [0.0];
# src = InitSources(pos, f0, ot, dt);
# # discretize spatial derivative operator
# vmax = maximum(v); vmin = minimum(v);
# fidMtx = InitFidMtx(nz, nx, ext, iflag, dz, dx, dt, vmax, vmin, f0, v);
# path = "/Users/wenlei/Desktop/acoustic.bin"
# MultiStepForward(path, src, fidMtx, tmax=1.0);
#
# using PyPlot, AcousticWave
# path = "/Users/wenlei/Desktop/acoustic.bin"
# pathout = "/Users/wenlei/Desktop/homo"
# waveAnim(path, pathout)
# sptPlot(path, 100)
# irx = collect(1:2:nx); irz = 5*ones(Int64,length(irx));
# pos_rec = hcat(irz, irx)
# shot = fullWfd2Shot(pos_rec, path)

ext: jl

max_stars_repo_name: WneleiGao/AcousticWaveLab

max_issues_repo_name: WneleiGao/AcousticWaveLab

max_forks_repo_name: WneleiGao/AcousticWaveLab


text:
module KumiKumiSlope

export part, base, color, entrances, exits, junctions, shape
export red, green, yellow, blue, orange, white
export phosphorescent_green, phosphorescent_yellow, clear_green, colorless
export set_100
export is_base, is_junction, is_leaf, is_root

include("common.jl")
include("parts.jl")
include("sets.jl")

end


text:
push!(LOAD_PATH, "../LF_common/")
using m_LF3d
include("../LF_common/sparse_LF3d.jl")
include("build_nabla2_x.jl")
include("build_nabla2_y.jl")
include("build_nabla2_z.jl")

const ⊗ = kron

function test_main( NN::Array{Int64} )

  println("NN = ", NN)

  Nx = NN[1]
  Ny = NN[2]
  Nz = NN[3]

  AA = [0.0, 0.0, 0.0]
  BB = [16.0, 16.0, 16.0]

  # Initialize LF
  LF = init_LF3d_c( NN, AA, BB, verbose=true )

  print("Building nabla2_x")
  @time Lx_v1 = build_nabla2_x(LF)

  print("Building nabla2_y")
  @time Ly_v1 = build_nabla2_y(LF)

  print("Building nabla2_z")
  @time Lz_v1 = build_nabla2_z(LF)

  #print("Forming Laplacian:")
  #@time ∇2 = Lx_v1 + Ly_v1

  #print("\n")

  #print("Using Kronecker product directly: nabla2_x")
  #@time Lx_v2 = LF.LFx.D2jl ⊗ speye(Ny) ⊗ speye(Nz)

  #print("Using Kronecker product directly: nabla2_y")
  #@time Ly_v2 = speye(Nx) ⊗ LF.LFy.D2jl ⊗ speye(Nz)

  #@printf("Using Kronecker product directly: nabla2_z")
  #@time Lz_v2 = speye(Nx) ⊗ speye(Ny) ⊗ LF.LFz.D2jl

  #print("Forming Laplacian:")
  #@time ∇2 = Lx_v2 + Ly_v2

end

@code_native test_main([3,3,3])

test_main( [80, 80, 80])


text:
function BodyCenterPointing(arg0::Frame, arg1::Ellipsoid)
    return BodyCenterPointing((Frame, Ellipsoid), arg0, arg1)
end

function get_attitude(obj::GroundPointing, arg0::FieldPVCoordinatesProvider, arg1::FieldAbsoluteDate, arg2::Frame)
    return jcall(obj, "getAttitude", FieldAttitude, (FieldPVCoordinatesProvider, FieldAbsoluteDate, Frame), arg0, arg1, arg2)
end

function get_attitude(obj::GroundPointing, arg0::PVCoordinatesProvider, arg1::AbsoluteDate, arg2::Frame)
    return jcall(obj, "getAttitude", Attitude, (PVCoordinatesProvider, AbsoluteDate, Frame), arg0, arg1, arg2)
end

function get_body_frame(obj::GroundPointing)
    return jcall(obj, "getBodyFrame", Frame, ())
end

function get_target_pv(obj::BodyCenterPointing, arg0::FieldPVCoordinatesProvider, arg1::FieldAbsoluteDate, arg2::Frame)
    return jcall(obj, "getTargetPV", TimeStampedFieldPVCoordinates, (FieldPVCoordinatesProvider, FieldAbsoluteDate, Frame), arg0, arg1, arg2)
end

function get_target_pv(obj::BodyCenterPointing, arg0::PVCoordinatesProvider, arg1::AbsoluteDate, arg2::Frame)
    return jcall(obj, "getTargetPV", TimeStampedPVCoordinates, (PVCoordinatesProvider, AbsoluteDate, Frame), arg0, arg1, arg2)
end


text:
#=
PSD cone
cuts are V*V^T for eigenvectors V corresponding to negative eigenvalues of matrix point X
=#
function get_init_cuts(cone::MOI.PositiveSemidefiniteConeTriangle)
    cuts = Vector{Float64}[]
    L = cone.side_dimension
    dim = MOI.dimension(cone)

    # diagonal variables are nonnegative
    for i in 1:L
        cut = zeros(dim)
        cut[mat_to_vec_idx(i, i)] = 1.0
        push!(cuts, cut)
    end

    # 2x2 principal minors satisfy linearizations of PSD/RSOC condition
    for i in 1:L, j in 1:(i - 1)
        cut = zeros(dim)
        cut[mat_to_vec_idx(i, i)] = 1.0
        cut[mat_to_vec_idx(j, j)] = 1.0
        ij_idx = mat_to_vec_idx(i, j)
        cut[ij_idx] = 2.0
        push!(cuts, cut)
        cut = copy(cut)
        cut[ij_idx] = -2.0
        push!(cuts, cut)
    end

    return cuts
end

function get_sep_cuts(
    x::Vector{Float64},
    cone::MOI.PositiveSemidefiniteConeTriangle,
    tol::Float64,
)
    cuts = Vector{Float64}[]
    L = cone.side_dimension
    X = Matrix{Float64}(undef, L, L)
    vec_to_mat_noscale!(X, x)

    # TODO try only getting negative eigenvalues (what lower bound?) vs smallest k eigenvalues
    F = eigen!(Symmetric(X, :L), 1:min(L, 5))

    for k in eachindex(F.values)
        if F.values[k] >= -0.1 * tol
            continue
        end

        V = F.vectors[:, k]
        W = V * V'
        cut = similar(x)
        mat_to_vec_scale!(cut, W)

        if dot(x, cut) < -tol
            push!(cuts, cut)
        end
    end

    return cuts
end


text:
d = [ 0 1 2 3 4; 5 0 6 7 8; 9 10 0 11 12; 13 14 15 0 16; 17 18 19 20 0]
z = floyd_warshall_shortest_paths(g3, d)
@test z.dists[3,:][:] == [7, 6, 0, 11, 27]
@test z.parents[3,:][:] == [2, 3, 0, 3, 4]

@test enumerate_paths(z)[2][2] == []
@test enumerate_paths(z)[2][4] == enumerate_paths(z,2)[4] == enumerate_paths(z,2,4) == [2,3,4]


text:
#### This example is a basic test of the new variables and factors
#### that are added in this repo. The example is derived from
#### the hexagonal test example.

using Caesar, RoME
using Caesar_VariableFactorExamples # Your new variable/factor repository
# Using plotting for experiment validation
using RoMEPlotting

# 1. Init factor graph
fg = initfg()

# 2. Add variables
addNode!(fg, :x0, ExamplePose2)
addNode!(fg, :x1, ExamplePose2)
addNode!(fg, :x2, ExamplePose2)
addNode!(fg, :x3, ExamplePose2)

# 3. Add factors
# 3a. Add a new test prior
# 3b. Add new types of odometry factors.

# 4. Solve graph
batchSolve!(fg)

# 5. Graph solution - assuming that you have this open in Atom.
drawPoses(fg)


text:
export athenasetup, athenaquery
"""
  athenasetup(bucket)

Create AWS Athena database and table for SQL queries.

# Arguments
- `bucket::String`: Name of AWS bucket to store query results
"""
function athenasetup(aws::AWSConfig,bucket::String;
                    database::String="scedcindex",
                    table::String="scedc_parquet",
                    clean::Bool=true)

    # make sure a bucket exists for queries
    if bucket ∉ s3_list_buckets(aws)
        s3_create_bucket(aws,bucket)
    end

    # create database
    println("Creating database '$database' ",now())
    Athena.start_query_execution("create database $database",
        Dict("ResultConfiguration"=>Dict("OutputLocation" => "s3://$bucket/queries/"),
        "ClientRequestToken"=>randstring(32)))

    # string to create table
    tablestr = """CREATE EXTERNAL TABLE IF NOT EXISTS $table ( `ms_filename` string,""" *
               """ `net` string, `sta` string, `seedchan` string, `location` string,""" *
               """ `lat` double, `lon` double, `sample_rate` double ) PARTITIONED BY""" *
               """ ( year int, year_doy string ) STORED AS PARQUET LOCATION""" *
               """ 's3://scedc-pds/continuous_waveforms/index/parquet/'""" *
               """ TBLPROPERTIES ("parquet.compress"="SNAPPY");"""

    sleep(1)

    # create table
    println("Creating table '$table' ",now())
    Athena.start_query_execution(tablestr,
       Dict("ResultConfiguration"=>Dict("OutputLocation" => "s3://$bucket/queries/"),
       "ClientRequestToken"=>randstring(32),
       "QueryExecutionContext" => Dict("Database" => database))
    )

    sleep(1)

    # update partitions
    println("Repairing table '$table' ",now())
    output = Athena.start_query_execution("MSCK REPAIR TABLE $table",
           Dict("ResultConfiguration"=>Dict("OutputLocation" => "s3://$bucket/queries/"),
           "ClientRequestToken"=>randstring(32),
           "QueryExecutionContext"=>Dict("Database" => database))
    )

    # wait until done
    finished = false
    while !finished
        queryresult = Athena.get_query_execution(output["QueryExecutionId"])
        if queryresult["QueryExecutionDetail"]["Status"]["State"] == "SUCCEEDED"
            finished = true
            println("REPAIR table $table finished ",now())
        elseif queryresult["QueryExecutionDetail"]["Status"]["State"] == "FAILED"
            error("REPAIR TABLE '$table' FAILED ",now())
        else
            sleep(10)
        end
    end
   return nothing
end
athenasetup(a...;b...) = athenasetup(global_aws_config(region="us-west-2"),a...;b...)
athenasetup(d::Dict,a...;b...) = athenasetup(global_aws_config(region=d[:region]),a...;b...)

"""
  athenaquery(bucket,query)

Use AWS Athena to query SCEDC-pds database.

Possible queries use these parameters:
- `year::String`: Year data was recorded, in YYYY format.
- `year_doy::String`: Year and day data was recorded, in YYYY_DDD format.
- `ms_filename::String`: Name of miniSEED file, e.g. SBCPSLOBHZ01_2020153.ms
- `net::String`: Network identifier, e.g. CI, AZ, etc..
- `seedchan::String`: SEED channel, e.g. BHZ, HHN, LHE, etc..
- `location::String`: Station location, e.g. 00, 01, 00, etc..
- `lat::double`: Station latitude, e.g. 34.5
- `lon::double`: Station longitude, e.g. -117
- `sample_rate::double`: Station sampling rate, in Hz

Logical operators available for queries include:
    - `Logical Operators`: `and`, `or`, `not`
    - `Comparison Operators`: `<`, `>`, `<=`, `>=`, `=`, `!=`
    - `Range Operator`: `between`
See https://docs.aws.amazon.com/athena/latest/ug/presto-functions.html for possible operators.

# Arguments
- `bucket::String`: Name of AWS bucket to store query results
- `query::String`: SQL query for SCEDC-pds database.

# Returns
Array of filepaths matching query in scedc-pds bucket.

# Example queries
1. Query all "HH[E,N,Z]" channels
query = "seedchan LIKE 'HH%'"

2. Query data from time range
query = "year_doy > '2016_150' and year_doy<'2016_157'"

3. Query lat/lon boundaries
query =  "lat between 34.00 and 34.50 and lon between -117 and -116"

4. Compound queries
query = "year = 2020 and seedchan LIKE 'BH%'"

"""
function athenaquery(aws::AWSConfig,bucket::String, query::String;
                    database::String="scedcindex",
                    table::String="scedc_parquet",
                    clean::Bool=true)

    # make sure a bucket exists for queries
    if bucket ∉ s3_list_buckets(aws)
        s3_create_bucket(aws,bucket)
    end

    # query the table
    output = Athena.start_query_execution(
        "select ms_filename from $table where $query order by ms_filename ",
        Dict("ResultConfiguration"=>Dict("OutputLocation" => "s3://$bucket/queries/"),
        "ClientRequestToken"=>randstring(32),
        "QueryExecutionContext" => Dict("Database" => database))
    )

    # check when query is finished
    finished = false
    while !finished
         queryresult = Athena.get_query_execution(output["QueryExecutionId"])
         if queryresult["QueryExecutionDetail"]["Status"]["State"] == "SUCCEEDED"
             finished = true
             println("QUERY '$query' finished ",now())
         elseif queryresult["QueryExecutionDetail"]["Status"]["State"] == "FAILED"
             error("QUERY '$query' FAILED ",now())
         else
             sleep(1)
         end
     end

    # read query results
    querypath = "queries/"*output["QueryExecutionId"] * ".csv"
    tmpfile = tempname()
    s3_get_file(aws,bucket,querypath,tmpfile)
    filelist = readdlm(tmpfile,String)
    filelist = filelist[2:end,1]
    if length(filelist) != 0
        println("$(length(filelist)) matching QUERIES ",now())
    else
        println("No matching QUERIES ",now())
    end

    if clean
        s3_delete(aws,bucket,querypath)
        s3_delete(aws,bucket,querypath*".metadata")
    end
    return scedcpath.(filelist)
end
athenaquery(a...;b...) = athenaquery(global_aws_config(region="us-west-2"),a...;b...)
athenaquery(d::Dict,a...;b...) = athenaquery(global_aws_config(region=d[:region]),a...;b...)


text:
# module MySQL_Octo.SQL

type MySQLAdapter <: QueryBuilders.AbstractAdapter
    handle::Union{Void, MySQLHandle}
end


text:
module test_02_OneScalarSignalWithUnit

using ModiaResult
using Unitful
using OrderedCollections
ModiaResult.@usingModiaPlot

t = range(0.0, stop=10.0, length=100)

result = OrderedDict{String,Any}()

result["time"] = t*u"s"
result["phi"]  = sin.(t)*u"rad"

println("\n... test_02_OneScalarSignalWithUnit.jl:\n")
ModiaResult.printResultInfo(result)

plot(result, "phi", heading="Sine(time)")

end

max_stars_repo_name: ModiaSim/ModiaResult

max_issues_repo_name: ModiaSim/ModiaResult

max_forks_repo_name: ModiaSim/ModiaResult


text:
# This file contains conversion functions for converting simplices with data to Scatter,
# LineSegments, and Mesh.

# This code is generated to make types specific enough.
# Using Union{Scatter, LineSegments, Mesh} does not work.
for T in (Scatter, LineSegments, Mesh)
    @eval begin
        function AbstractPlotting.convert_arguments(
            ::Type{T}, sx::AbstractSimplex, data::AbstractVector
        ) where {T<:$T}
            return convert_arguments(T, [sx], data)
        end

        function AbstractPlotting.convert_arguments(
            ::Type{T}, elem::AbstractChainElement, data::AbstractVector
        ) where {T<:$T}
            return convert_arguments(T, [simplex(elem)], data)
        end

        function AbstractPlotting.convert_arguments(
            ::Type{T}, chain::AbstractVector{<:AbstractChainElement}, data::AbstractVector
        ) where {T<:$T}
            return convert_arguments(T, simplex.(chain), data)
        end
    end
end

function AbstractPlotting.convert_arguments(
    ::Type{T}, chain::AbstractVector{<:AbstractSimplex}, data::AbstractVector
) where {T<:Scatter}
    _data, = convert_arguments(PointBased(), data)
    points = unique!(collect(Iterators.flatten(_data[sx] for sx in chain)))
    return convert_arguments(T, points)
end
function AbstractPlotting.convert_arguments(
    ::Type{T}, chain::AbstractVector{<:AbstractSimplex}, data::AbstractVector
) where {T<:LineSegments}
    _data, = convert_arguments(PointBased(), data)
    segs = NTuple{2,eltype(data)}[]
    for sx in chain
        for vs in IterTools.subsets(vertices(sx), Val(2))
            push!(segs, getindex.(Ref(data), vs))
        end
    end
    unique!(segs)
    return convert_arguments(T, collect(Iterators.flatten(segs)))
end
function AbstractPlotting.convert_arguments(
    ::Type{T}, chain::AbstractVector{<:AbstractSimplex}, data::AbstractVector
) where {T<:Mesh}
    _data, = convert_arguments(PointBased(), data)
    faces = typeof(GeometryBasics.GLTriangleFace(1, 2, 3))[]
    for sx in chain
        for vs in IterTools.subsets(vertices(sx), Val(3))
            push!(faces, vs)
        end
    end
    return convert_arguments(T, GeometryBasics.Mesh(_data, faces))
end


text:
"""
    splitobs(n::Int; at) -> Tuple

Compute the indices for two or more disjoint subsets of
the range `1:n` with splits given by `at`.

# Examples

```julia
julia> splitobs(100, at=0.7)
(1:70, 71:100)

julia> splitobs(100, at=(0.1, 0.4))
(1:10, 11:50, 51:100)
```
"""
splitobs(n::Int; at) = _splitobs(n, at)

_splitobs(n::Int, at::Integer) = _splitobs(n::Int, at / n)
_splitobs(n::Int, at::NTuple{N, <:Integer}) where {N} = _splitobs(n::Int, at ./ n)

_splitobs(n::Int, at::Tuple{}) = (1:n,)

function _splitobs(n::Int, at::AbstractFloat)
    0 <= at <= 1 || throw(ArgumentError("the parameter \"at\" must be in interval (0, 1)"))
    n1 = clamp(round(Int, at*n), 0, n)
    (1:n1, n1+1:n)
end

function _splitobs(n::Int, at::NTuple{N,<:AbstractFloat}) where N
    at1 = first(at)
    a, b = _splitobs(n::Int, at1)
    n1 = a.stop
    n2 = b.stop
    at2 = Base.tail(at) .* n ./ (n2 - n1)
    rest = map(x -> n1 .+ x, _splitobs(n2-n1, at2))
    return (a, rest...)
end

"""
    splitobs(data; at, shuffle=false) -> Tuple

Split the `data` into multiple subsets proportional to the
value(s) of `at`.

If `shuffle=true`, randomly permute the observations before splitting.

Supports any datatype implementing the [`numobs`](@ref) and
[`getobs`](@ref) interfaces.

# Examples

```julia
# A 70%-30% split
train, test = splitobs(X, at=0.7)

# A 50%-30%-20% split
train, val, test = splitobs(X, at=(0.5, 0.3))

# A 70%-30% split with multiple arrays and shuffling
train, test = splitobs((X, y), at=0.7, shuffle=true)
Xtrain, Ytrain = train
```
"""
function splitobs(data; at, shuffle::Bool=false)
    if shuffle
        data = shuffleobs(data)
    end
    n = numobs(data)
    return map(idx -> obsview(data, idx), splitobs(n; at))
end

max_stars_repo_name: JuliaML/MLBase

max_issues_repo_name: JuliaML/MLBase

max_forks_repo_name: JuliaML/MLBase


text:
using LinearAlgebra;
using GenericSchur;

function diag_fun!(fun, A::Array{T, 1}, B::Array{T, 1}, C::Array{T, 2}) where T

    (m, n) = size(C)

    for j = 1 : n
        for i = 1 : m
            C[i,j] = fun(A[i], B[j], 0, 0) .* C[i,j];
        end
    end

    return C;
end

function diag_fun(fun, A::Array{T, 1}, B::Array{T, 1}, C::Array{T, 2}) where T
    F = copy(C);
    diag_fun!(fun, A, B, F)
end

function isdiagonal(A::Matrix)
    (m, n) = size(A)
    for i = 1 : m
        for j = 1 : n
            if i != j && A[i,j] != 0
                return false
            end
        end
    end

    return true
end

function diag_fun(fun, A::Matrix, B::Matrix, C::Matrix)
    F = similar(C);

    isdiagA = isdiagonal(A);
    isdiagB = isdiagonal(B);

    # Computes fun{A, B}(C) when at least one among A and B is diagonal
    if isdiagA && isdiagB
        F .= C;
        diag_fun!(fun, diag(A), diag(B), C)
        return F;
    end

    if isdiagA
        dA = diag(A);
        for j = 1:size(A, 1)
            # g = @(y, k) fun(dA(j), y, 0, k);
            g = (y,k) -> fun(dA[j], y, 0, k);
            # funB = @(y, k) arrayfun(g, y, k * ones(size(y)));
            funB = (y, k) -> map(g, y, k * ones(size(y)));
            F[j, :] = C[j, :] * funm(B, funB);
        end
        return F;
    end

    if isdiagB
        dB = diag(B);
        for j = 1:size(B, 1)
            # g = @(x, k) fun(x, dB(j), k, 0);
            g = (x, k) -> fun(x, dB(j), k, 0);
            # funA = @(x, k) arrayfun(g, x, k * ones(size(x)));
            funA = (x, k) -> map(g, x, k * ones(size(x)));
            F[:, j] = funm(A, funA) * C(:, j);
        end
        return F;
    end

    ta = Threads.@spawn eigen(A);
    tb = Threads.@spawn eigen(B);

    A, VA = fetch(ta)
    B, VB = fetch(tb)

    F = VA * diag_fun(fun, A, B, VA \ C * VB) / VB;
end

#----------------------------------- Auxiliary functions-----------------------------------

"""
    fun2m_atom(f, A::Matrix, B::Matrix, C::Matrix, treeA::BivMatTree, treeB::BivMatTree, meth::Algoroithm, min_digits::Int64, user_function = nothing)

    Bivariate function of triangular coefficients with nearly constant diagonals.
"""
function fun2m_atom(f, A::Matrix, B::Matrix, C::Matrix, treeA::BivMatTree, treeB::BivMatTree, meth::Algorithm, min_digits::Int64, user_function = nothing)
    digits = 0;

    nA = size(A, 1); nB = size(B, 1);

    if meth == Taylor
        F, digits = taylor_eval(f, A, B, C, eps());
    elseif meth == Diag
        uh = eps() / treeA.kV / treeB.kV;
        d_uh = convert(Int64, max(ceil(log10(1 / uh)), min_digits));

        F, digits = fun2m_atom_pad(treeA.V, treeA.D, treeB.V, treeB.D, C, f, d_uh);
    elseif meth == DiagNoHp # diagonalization without high precision
        VA = treeA.V; DA = treeA.D;
        VB = treeB.V; DB = treeB.D;
        F = VA \ C * VB;
        diag_fun!(f, DA, DB, F);
        F = VA * F / VB;
    elseif meth == User # handle function provided by the user
        F = user_function(A, B, C);
    end

    return F, digits;
end

"""
    fun2m_atom_pad!(VA, DA, VB, DB, C, fun, d_uh)

    Multiprecision perturb-and-diagonalize for evaluating bivariate
    functions of matrices f{A, B^T}(C) for upper triangular A, B; it
    is expected that the eigenvector basis and the eigenvalues are
    precomputed in VA, DA, VB, DB.
"""
function fun2m_atom_pad(VA, DA, VB, DB, C, fun, d_uh)
    F = copy(C)
    F, d_uh = fun2m_atom_pad!(VA, DA, VB, DB, F, fun, d_uh)
end

"""
    fun2m_atom_pad!(VA, DA, VB, DB, C, fun, d_uh)

    Multiprecision perturb-and-diagonalize for evaluating bivariate
    functions of matrices f{A, B^T}(C) for upper triangular A, B; it
    is expected that the eigenvector basis and the eigenvalues are
    precomputed in VA, DA, VB, DB.

    The evaluation is performed in-place, overwriting C.
"""
function fun2m_atom_pad!(VA, DA, VB, DB, C, fun, d_uh)

    if d_uh <= 18
        VA = UpperTriangular(convert(Matrix{ComplexF64}, VA))
        VB = UpperTriangular(convert(Matrix{ComplexF64}, VB))
        C  = convert(Matrix{ComplexF64}, C)
        DA = convert(Array{ComplexF64}, DA)
        DB = convert(Array{ComplexF64}, DB)

        C = VA \ C * VB;
        diag_fun!(fun, DA, DB, C)

        C = VA * C / VB
    else
        C = with_digits(d_uh) do
            VA = UpperTriangular(mp(VA, d_uh))
            VB = UpperTriangular(mp(VB, d_uh))
            DA = mp(DA, d_uh)
            DB = mp(DB, d_uh)
            C = mp(C, d_uh)

            C = VA \ C;
            # C = bf_mldivide!(VA, C)
            # C = C * VB;
            C = bf_matmulr!(C, VB)
            diag_fun!(fun, DA, DB, C)
            C = VA * C;
            C = C / VB

            convert(Matrix{ComplexF64}, C);
        end
    end

    return C, d_uh
end

function bf_mldivide!(VA, C)
    # In place C = VA \ C, with upper triangular VA, stored in C.
    n = size(VA, 1)
    m = size(C, 2)

    for i = n : -1 : 1
        for j = 1 : m
            for k = i + 1 : n
                C[i,j] = C[i,j] - VA[i,k] * C[k,j]
            end
            C[i,j] = C[i,j] / VA[i,i]
        end
    end

    C
end

function bf_matmulr!(A, VB)
    # In place matmul C = A*VB, with upper triangular VB, stored in A.
    m = size(A, 1);
    n = size(A, 2);

    for j = n : -1 : 1
        for k = 1 : m
            A[k, j] = A[k, j] * VB[j,j];
        end

        for i = j-1 : -1 : 1
            for k = 1 : m
                A[k, j] = A[k, j] + VB[i,j] * A[k, i];
            end
        end
    end

    A
end

# function bf_matmulr2!(A, VB)
#     # In place matmul C = A*VB, with upper triangular VB, stored in A.
#     m = size(A, 1);
#     n = size(A, 2);
#     p = size(VB, 2);

#     tmp = complex(MPFR_TLS.BigFloat());
#     tmp1 = MPFR_TLS.BigFloat();
#     tmp2 = MPFR_TLS.BigFloat();

#     r = MPFR_TLS.ROUNDING_MODE[];

#     for j = p : -1 : 1
#         for k = 1 : m
#             # A[k, j] = A[k, j] * VB[j,j];
#             ccall((:mpfr_mul, :libmpfr), Int32, (Ref{MPFR_TLS.BigFloat}, Ref{MPFR_TLS.BigFloat}, Ref{MPFR_TLS.BigFloat}, MPFR_TLS.MPFRRoundingMode),
#                 tmp1, real(VB[j,j]), real(A[k,j]) ,r);
#             ccall((:mpfr_mul, :libmpfr), Int32, (Ref{MPFR_TLS.BigFloat}, Ref{MPFR_TLS.BigFloat}, Ref{MPFR_TLS.BigFloat}, MPFR_TLS.MPFRRoundingMode),
#                 tmp2, imag(VB[j,j]), imag(A[k,j]) ,r);
#             tmp1 = real(A[k,j])
#             ccall((:mpfr_sub, :libmpfr), Int32, (Ref{MPFR_TLS.BigFloat}, Ref{MPFR_TLS.BigFloat}, Ref{MPFR_TLS.BigFloat}, MPFR_TLS.MPFRRoundingMode),
#                 real(A[k,j]), tmp1, tmp2, r);

#             ccall((:mpfr_mul, :libmpfr), Int32, (Ref{MPFR_TLS.BigFloat}, Ref{MPFR_TLS.BigFloat}, Ref{MPFR_TLS.BigFloat}, MPFR_TLS.MPFRRoundingMode),
#                 tmp1, real(VB[j,j]), imag(A[k,j]) ,r);
#             ccall((:mpfr_add, :libmpfr), Int32, (Ref{MPFR_TLS.BigFloat}, Ref{MPFR_TLS.BigFloat}, Ref{MPFR_TLS.BigFloat}, MPFR_TLS.MPFRRoundingMode),
#                 tmp2, imag(VB[j,j]), tmp1, r);
#             ccall((:mpfr_add, :libmpfr), Int32, (Ref{MPFR_TLS.BigFloat}, Ref{MPFR_TLS.BigFloat}, Ref{MPFR_TLS.BigFloat}, MPFR_TLS.MPFRRoundingMode),
#                 imag(A[k,j]), tmp1, tmp2, r);
#         end

#         for i = j-1 : -1 : 1
#             for k = 1 : m
#                 # A[k, j] = A[k, j] + VB[i,j] * A[k, i];
#                 ccall((:mpfr_mul, :libmpfr), Int32, (Ref{MPFR_TLS.BigFloat}, Ref{MPFR_TLS.BigFloat}, Ref{MPFR_TLS.BigFloat}, MPFR_TLS.MPFRRoundingMode),
#                     tmp1, real(VB[i,j]), real(A[k,i]), r);
#                 ccall((:mpfr_mul, :libmpfr), Int32, (Ref{MPFR_TLS.BigFloat}, Ref{MPFR_TLS.BigFloat}, Ref{MPFR_TLS.BigFloat}, MPFR_TLS.MPFRRoundingMode),
#                     tmp2, imag(VB[i,j]), imag(A[k,i]), r);
#                 ccall((:mpfr_sub, :libmpfr), Int32, (Ref{MPFR_TLS.BigFloat}, Ref{MPFR_TLS.BigFloat}, Ref{MPFR_TLS.BigFloat}, MPFR_TLS.MPFRRoundingMode),
#                     real(tmp), tmp1, tmp2, r);

#                 ccall((:mpfr_mul, :libmpfr), Int32, (Ref{MPFR_TLS.BigFloat}, Ref{MPFR_TLS.BigFloat}, Ref{MPFR_TLS.BigFloat}, MPFR_TLS.MPFRRoundingMode),
#                     tmp1, real(VB[i,j]), imag(A[k,i]), r);
#                 ccall((:mpfr_mul, :libmpfr), Int32, (Ref{MPFR_TLS.BigFloat}, Ref{MPFR_TLS.BigFloat}, Ref{MPFR_TLS.BigFloat}, MPFR_TLS.MPFRRoundingMode),
#                     tmp2, imag(VB[i,j]), real(A[k,i]), r);
#                 ccall((:mpfr_add, :libmpfr), Int32, (Ref{MPFR_TLS.BigFloat}, Ref{MPFR_TLS.BigFloat}, Ref{MPFR_TLS.BigFloat}, MPFR_TLS.MPFRRoundingMode),
#                     imag(tmp), tmp1, tmp2, r);

#                 ccall((:mpfr_add, :libmpfr), Int32, (Ref{MPFR_TLS.BigFloat}, Ref{MPFR_TLS.BigFloat}, Ref{MPFR_TLS.BigFloat}, MPFR_TLS.MPFRRoundingMode),
#                     real(A[k,j]), real(tmp), real(A[k,j]) ,r);
#                 ccall((:mpfr_add, :libmpfr), Int32, (Ref{MPFR_TLS.BigFloat}, Ref{MPFR_TLS.BigFloat}, Ref{MPFR_TLS.BigFloat}, MPFR_TLS.MPFRRoundingMode),
#                     imag(A[k,j]), imag(tmp), imag(A[k,j]) ,r);
#             end
#         end
#     end

#     A
# end

function fun2m_taylor(fun, x, y)
    # FUN2M_TAYLOR

    rx = maximum(abs(x));
    ry = maximum(abs(y));

    # Start by computing a degree 2 Taylor expansion
    C = zeros(3);

    for i = 1 : 3
        for j = 1 : 3
            C[i,j] = fun(0, 0, i-1, j-1) / factorial(i-1) / factorial(j-1);
        end
    end

    can_stop = false;

    dx = 2;
    dy = 2;

    tx = rx^dx;
    ty = ry^dy;

    vx = [ 1, rx, rx.^2 ];
    vy = [ 1, ry, ry.^2 ];

    fx = 1 / 2;
    fy = 1 / 2;

    tol = eps();

    while !can_stop
        can_stop = true;

        res_x = norm( (C[end, :] .* vy) * tx, Inf);
        res_y = norm( (C[:, end] .* vx) * ty, Inf);

        if res_x > tol * maximum(abs(C(:)))
            dx = dx + 1;
            C = [ C ; zeros(1, size(C, 2)) ];

            fx = fx / dx;

            for i = 1 : size(C, 2)
                C[end, i] = fun(0, 0, size(C,1)-1, i-1) * fx / factorial(i-1);
            end

            tx = tx * rx;
            vx = [ vx  vx[end] * rx ];
            can_stop = false;
        end

        if res_y > tol * maximum(abs(C))
            dy = dy + 1;
            C = [ C zeros(size(C, 1), 1) ];

            fy = fy / dy;

            for i = 1 : size(C, 1)
                C[i, end] = fun(0, 0, i-1, size(C,2)-1) * fy / factorial(i-1);
            end

            ty = ty * ry;
            vy = [ vy vy[end]*ry ];
            can_stop = false;
        end

        if can_stop
            # FIXME: We may want to check the maximum of the derivative over
            # the rectangle [-rx, rx] x [-ry, ry], to see if we can really stop
            C = C[1:end-1, 1:end-1];
        end

        return C
    end

end

max_stars_repo_name: numpi/BivMatFun

max_issues_repo_name: numpi/BivMatFun

max_forks_repo_name: numpi/BivMatFun


text:
"""
    attribute_value(x)

This method may be implemented to specify a printed representation
suitable for use within a quoted attribute value. By default, the print
representation of an object is used, and then propertly escaped. There
are a few overrides that we provide.

* The elements of a `Tuple` or `AbstractArray` object are printed,
  with a space between each item.

* The `Pair`, `NamedTuple`, and `Dict` objects are treated as if
  they are CSS style elements, with a colon between key and value,
  each pair delimited by a semi-colon.

* The `Bool` object, which has special treatment for bare inside_tag,
  is an error when used within a quoted attribute.

If an object is wrapped with `HTML` then it is included in the quoted
attribute value as-is, without inspection or escaping.
"""
attribute_value(x::AbstractString) = x
attribute_value(x::Number) = x
attribute_value(x::Symbol) = x
attribute_value(x::Nothing) = ""
attribute_value(x::Bool) =
  throw("Boolean used within a quoted attribute.")

function attribute_value(xs::Union{Tuple, AbstractArray, Base.Generator})
    Reprint() do io::IO
        prior = false
        for x in xs
            if prior
                print(io, " ")
            end
            print(io, attribute_value(x))
            prior = true
        end
    end
end

function attribute_values(xs)
    Reprint() do io::IO
        prior = false
        for (key, value) in xs
            name = normalize_attribute_name(key)
            if prior
                print(io, "; ")
            end
            print(io, name)
            print(io, ": ")
            print(attribute_value(value))
            prior = true
        end
        print(io, ";")
    end
end

attribute_value(pair::Pair) = attribute_values((pair,))
attribute_value(items::Dict) = attribute_values(items)
attribute_value(items::NamedTuple) = attribute_values(pairs(items))
attribute_value(items::Tuple{Pair, Vararg{Pair}}) = attribute_values(items)

"""
    content(x)

This method may be implemented to specify a printed representation
suitable for `text/html` output. As a special case, if the result is
wrapped with `HTML`, then it is passed along as-is. Otherwise, the
`print` representation of the resulting value is escaped. By default
`AbstractString`, `Number` and `Symbol` values are printed and escaped.
The elements of `Tuple` and `AbstractArray` are concatinated and then
escaped. If a method is not implemented for a given object, then we
attempt to `show` it via `MIME"text/html"`.
"""
content(x) = Render(x)
content(x::AbstractString) = x
content(x::Number) = x
content(x::Symbol) = x
content(x::Nothing) = ""
content(xs...) = content(xs)

function content(xs::Union{Tuple, AbstractArray, Base.Generator})
    Reprint() do io::IO
        for x in xs
            print(io, content(x))
        end
    end
end

#-------------------------------------------------------------------------
"""
    attribute_pair(name, value)

Wrap and escape attribute name and pair within a single-quoted context
so that it is `showable("text/html")`. It's assumed that the attribute
name has already been normalized.

If an attribute value is `Bool` or `Nothing`, then special treatment is
provided. If the value is `false` or `nothing` then the entire pair is
not printed.  If the value is `true` than an empty string is produced.
"""

no_content = Reprint(io::IO -> nothing)

function attribute_pair(name, value)
    Reprint() do io::IO
        print(io, " ")
        print(io, name)
        print(io, Bypass("='"))
        print(io, attribute_value(value))
        print(io, Bypass("'"))
    end
end

function attribute_pair(name, value::Bool)
    if value == false
        return no_content
    end
    Reprint() do io::IO
        print(io, " ")
        print(io, name)
        print(io, Bypass("=''"))
    end
end

attribute_pair(name, value::Nothing) = no_content

"""
    inside_tag(value)

Convert Julian object into a serialization of attribute pairs,
`showable` via `MIME"text/html"`. The default implementation of this
delegates value construction of each pair to `attribute_pair()`.
"""
function inside_tag(value::Pair)
    name = normalize_attribute_name(value.first)
    return attribute_pair(name, value.second)
end

function inside_tag(xs)
    Reprint() do io::IO
        for (key, value) in xs
            name = normalize_attribute_name(key)
            print(io, attribute_pair(name, value))
        end
    end
end

inside_tag(values::NamedTuple) =
    inside_tag(pairs(values))

"""
    rawtext(context, value)

Wrap a string value that occurs with RAWTEXT, SCRIPT and other element
context so that it is `showable("text/html")`. The default
implementation ensures that the given value doesn't contain substrings
illegal for the given context.
"""
function rawtext(context::Symbol, value::AbstractString)
    if occursin("</$context>", lowercase(value))
        throw(DomainError(repr(value), "  Content of <$context> cannot " *
            "contain the end tag (`</$context>`)."))
    end
    if context == :script && occursin("<!--", value)
        # this could be slightly more nuanced
        throw(DomainError(repr(value), "  Content of <$context> should " *
            "not contain a comment block (`<!--`) "))
    end
    return Bypass(value)
end


text:
module DiffGPFlux

    using Reexport
    @reexport using ForwardDiff
    @reexport using Distributions
    @reexport using LinearAlgebra

    include("exact_gp.jl")
    include("kernel.jl")
    include("utils.jl")
    include("mean.jl")
    include("optimizers.jl")

    export GaussProcess, negloglik, predict, ∇ₓpredict
    export rbf, ∇ₓrbf, linear, ∇ₓlinear, mixed_effect_kernel
    export ConstMean, LinearMean
    export seed_duals, update!, grads
    export Descent, Momentum, Nesterov, RMSProp, ADAM, ADAGrad

end # module

max_stars_repo_name: HamletWantToCode/DiffGPFlux

max_issues_repo_name: HamletWantToCode/DiffGPFlux


text:
@testset "765.couples-holding-hands.jl" begin
    @test min_swaps_couples([3, 2, 0, 1]) == 0
    @test min_swaps_couples([0, 2, 1, 3]) == 1
    @test min_swaps_couples([24, 58, 43, 7, 17, 3, 6, 54, 27, 14, 18, 57, 53, 50, 12, 42,
                             49, 1, 2, 39, 47, 56, 38, 8, 45, 22, 36, 0, 30, 52, 23, 40, 16,
                             4, 35, 55, 10, 34, 11, 48, 28, 46, 32, 15, 29, 13, 9, 59, 44,
                             51, 5, 21, 33, 26, 25, 41, 19, 31, 37, 20]) == 28
end


text:
type CurrentPlot
    nullableplot::Nullable{AbstractPlot}
end
const CURRENT_PLOT = CurrentPlot(Nullable{AbstractPlot}())

isplotnull() = isnull(CURRENT_PLOT.nullableplot)

function current()
    if isplotnull()
        error("No current plot/subplot")
    end
    get(CURRENT_PLOT.nullableplot)
end
current(plot::AbstractPlot) = (CURRENT_PLOT.nullableplot = Nullable(plot))

# ---------------------------------------------------------

Base.string(plt::Plot) = "Plot{$(plt.backend) n=$(plt.n)}"
Base.print(io::IO, plt::Plot) = print(io, string(plt))
Base.show(io::IO, plt::Plot) = print(io, string(plt))

getplot(plt::Plot) = plt
getplotargs(plt::Plot, idx::Int = 1) = plt.plotargs
convertSeriesIndex(plt::Plot, n::Int) = n

# ---------------------------------------------------------

"""
The main plot command.  Use `plot` to create a new plot object, and `plot!` to add to an existing one:

```
    plot(args...; kw...)                  # creates a new plot window, and sets it to be the current
    plot!(args...; kw...)                 # adds to the `current`
    plot!(plotobj, args...; kw...)        # adds to the plot `plotobj`
```

There are lots of ways to pass in data, and lots of keyword arguments... just try it and it will likely work as expected.
When you pass in matrices, it splits by columns.  See the documentation for more info.
"""

# this creates a new plot with args/kw and sets it to be the current plot
function plot(args...; kw...)
    pkg = backend()
    d = KW(kw)
    preprocessArgs!(d)
    dumpdict(d, "After plot preprocessing")

    plotargs = merge(d, getPlotArgs(pkg, d, 1))
    dumpdict(plotargs, "Plot args")
    plt = _create_plot(pkg, plotargs)  # create a new, blank plot

    delete!(d, :background_color)
    plot!(plt, args...; kw...)  # add to it
end

# this adds to the current plot, or creates a new plot if none are current
function  plot!(args...; kw...)
    local plt
    try
        plt = current()
    catch
        return plot(args...; kw...)
    end
    plot!(current(), args...; kw...)
end

# this adds to a specific plot... most plot commands will flow through here
function plot!(plt::Plot, args...; kw...)
    d = KW(kw)
    userkw = preprocessArgs!(d)

    # for plotting recipes, swap out the args and update the parameter dictionary
    args = RecipesBase.apply_recipe(d, userkw, args...)
    _add_markershape(d)

    dumpdict(d, "After plot! preprocessing")
    warnOnUnsupportedArgs(plt.backend, d)

    # just in case the backend needs to set up the plot (make it current or something)
    _before_add_series(plt)

    # # grouping
    groupby = if haskey(d, :group)
        extractGroupArgs(d[:group], args...)
    else
        nothing
    end

    # merge plot args
    if !haskey(d, :subplot)
        for k in keys(_plotDefaults)
            if haskey(d, k)
                plt.plotargs[k] = d[k]
            end
        end
        # merge!(plt.plotargs, d)
        handlePlotColors(plt.backend, plt.plotargs)
    end

    _add_series(plt, d, groupby, args...)
    _add_annotations(plt, d)

    warnOnUnsupportedScales(plt.backend, d)

    # add title, axis labels, ticks, etc
    if !haskey(d, :subplot)
        merge!(plt.plotargs, d)
        # handlePlotColors(plt.backend, plt.plotargs)
        dumpdict(plt.plotargs, "Updating plot items")
        _update_plot(plt, plt.plotargs)
    end

    _update_plot_pos_size(plt, d)

    current(plt)

    # note: lets ignore the show param and effectively use the semicolon at the end of the REPL statement
    # # do we want to show it?
    if haskey(d, :show) && d[:show]
        gui()
    end

    plt
end

# handle the grouping
function _add_series(plt::Plot, d::KW, groupby::GroupBy, args...)
    starting_n = plt.n
    for (i, glab) in enumerate(groupby.groupLabels)
        tmpd = copy(d)
        tmpd[:numUncounted] = plt.n - starting_n
        _add_series(plt, tmpd, nothing, args...;
                    idxfilter = groupby.groupIds[i],
                    grouplabel = string(glab))
    end
end

filter_data(v::AVec, idxfilter::AVec{Int}) = v[idxfilter]
filter_data(v, idxfilter) = v

function filter_data!(d::KW, idxfilter)
    for s in (:x, :y, :z)
        d[s] = filter_data(get(d, s, nothing), idxfilter)
    end
end

function _replace_linewidth(d::KW)
    # get a good default linewidth... 0 for surface and heatmaps
    if get(d, :linewidth, :auto) == :auto
        d[:linewidth] = (get(d, :linetype, :path) in (:surface,:heatmap,:image) ? 0 : 1)
    end
end

# no grouping
function _add_series(plt::Plot, d::KW, ::Void, args...;
                     idxfilter = nothing,
                     grouplabel = "")

    # get the list of dictionaries, one per series
    dumpdict(d, "before process_inputs")
    process_inputs(plt, d, args...)
    dumpdict(d, "after process_inputs")

    if idxfilter != nothing
        # add the group name as the label if there isn't one passed in
        get!(d, :label, grouplabel)
        # filter the data
        filter_data!(d, idxfilter)
    end
    # dumpdict(d,"",true)

    seriesArgList, xmeta, ymeta = build_series_args(plt, d) #, idxfilter)
    # seriesArgList, xmeta, ymeta = build_series_args(plt, groupargs..., args...; d...)

    # if we were able to extract guide information from the series inputs, then update the plot
    # @show xmeta, ymeta
    updateDictWithMeta(d, plt.plotargs, xmeta, true)
    updateDictWithMeta(d, plt.plotargs, ymeta, false)

    # now we can plot the series
    for (i,di) in enumerate(seriesArgList)
        plt.n += 1

        if !stringsSupported() && di[:linetype] != :pie
            setTicksFromStringVector(plt, d, di, "x")
            setTicksFromStringVector(plt, d, di, "y")
            setTicksFromStringVector(plt, d, di, "z")
        end

        # remove plot args
        for k in keys(_plotDefaults)
            delete!(di, k)
        end

        # merge in plotarg_overrides
        plotarg_overrides = pop!(di, :plotarg_overrides, nothing)
        if plotarg_overrides != nothing
            merge!(plt.plotargs, plotarg_overrides)
        end
        # dumpdict(plt.plotargs, "pargs", true)

        dumpdict(di, "Series $i")

        _replace_linewidth(di)

        _add_series(plt.backend, plt, di)
    end
end

# --------------------------------------------------------------------

function get_indices(orig, labels)
    Int[findnext(labels, l, 1) for l in orig]
end

function setTicksFromStringVector(plt::Plot, d::KW, di::KW, letter)
    sym = symbol(letter)
    ticksym = symbol(letter * "ticks")
    pargs = plt.plotargs
    v = di[sym]

    # do we really want to do this?
    typeof(v) <: AbstractArray || return
    isempty(v) && return
    trueOrAllTrue(_ -> typeof(_) <: AbstractString, v) || return

    # compute the ticks and labels
    ticks, labels = if ticksType(pargs[ticksym]) == :ticks_and_labels
        # extend the existing ticks and labels. only add to labels if they're new!
        ticks, labels = pargs[ticksym]
        newlabels = filter(_ -> !(_ in labels), unique(v))
        newticks = if isempty(ticks)
            collect(1:length(newlabels))
        else
            maximum(ticks) + collect(1:length(newlabels))
        end
        ticks = vcat(ticks, newticks)
        labels = vcat(labels, newlabels)
        ticks, labels
    else
        # create new ticks and labels
        newlabels = unique(v)
        collect(1:length(newlabels)), newlabels
    end

    d[ticksym] = ticks, labels
    plt.plotargs[ticksym] = ticks, labels

    # add an origsym field so that later on we can re-compute the x vector if ticks change
    origsym = symbol(letter * "orig")
    di[origsym] = v
    di[sym] = get_indices(v, labels)

    # loop through existing plt.seriesargs and adjust indices if there is an origsym key
    for sargs in plt.seriesargs
        if haskey(sargs, origsym)
            # TODO: might need to call the setindex function instead to trigger a plot update for some backends??
            sargs[sym] = get_indices(sargs[origsym], labels)
        end
    end
end

# --------------------------------------------------------------------

_before_add_series(plt::Plot) = nothing

# --------------------------------------------------------------------

# should we update the x/y label given the meta info during input slicing?
function updateDictWithMeta(d::KW, plotargs::KW, meta::Symbol, isx::Bool)
    lsym = isx ? :xlabel : :ylabel
    if plotargs[lsym] == default(lsym)
        d[lsym] = string(meta)
    end
end
updateDictWithMeta(d::KW, plotargs::KW, meta, isx::Bool) = nothing

# --------------------------------------------------------------------

annotations(::@compat(Void)) = []
annotations{X,Y,V}(v::AVec{@compat(Tuple{X,Y,V})}) = v
annotations{X,Y,V}(t::@compat(Tuple{X,Y,V})) = [t]
annotations(v::AVec{PlotText}) = v
annotations(v::AVec) = map(PlotText, v)
annotations(anns) = error("Expecting a tuple (or vector of tuples) for annotations: ",
                       "(x, y, annotation)\n    got: $(typeof(anns))")

function _add_annotations(plt::Plot, d::KW)
    anns = annotations(get(d, :annotation, nothing))
    if !isempty(anns)

        # if we just have a list of PlotText objects, then create (x,y,text) tuples
        if typeof(anns) <: AVec{PlotText}
            x, y = plt[plt.n]
            anns = Tuple{Float64,Float64,PlotText}[(x[i], y[i], t) for (i,t) in enumerate(anns)]
        end

        _add_annotations(plt, anns)
    end
end

# --------------------------------------------------------------------

function Base.copy(plt::Plot)
    backend(plt.backend)
    plt2 = plot(; plt.plotargs...)
    for sargs in plt.seriesargs
        sargs = filter((k,v) -> haskey(_seriesDefaults,k), sargs)
        plot!(plt2; sargs...)
    end
    plt2
end

# --------------------------------------------------------------------


text:
using Test
using Viznet

@testset "unit disk" begin
    @test rand_unitdisk(10, 1.8) isa UnitDisk
    ud = UnitDisk([(0.1, 0.02), (0.12, 0.03), (0.3, 0.3)], 0.1)
    @test length(ud) == 3
    @test (bonds(ud) |> length) == 1
    @test Base.show(stdout, MIME("text/html"), ud) == false
end

include("parallelogram.jl")
include("squarelattice.jl")

include("chimera.jl")


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

# TODO: function generate_output_atomic(f::Function, outputs, inputs)


text:
function DisparateImpactRemover(X, protected_attribute, select_grp)
	grps = levels(X[:,protected_attribute])
	for grp in grps
		select_grp = X[:,protected_attribute] .== grp
		category = X[select_grp,:]
		all_grps[grp] = nquantile(category[:,repair_column],10000)
	end
	medians = Dict()
	for index in 1:10000
		medians[index] = []
	end
	for index in 1:10000
		for grp in grps
			append!(medians[index],all_grps[grp][index])
		end
	end
	for grp in grps
		select_grp = X[:,protected_attribute] .== grp
		category = X[select_grp,:]
		for (index,var) in enumerate(category[:,select_grp])
			try
				category[index,select_grp]=sort!(medians[convert(Int64,floor(cdf(var)*10000))])[4]
            catch
				println("Couldn't replace")
			end
		end
	end
end

struct DisparateImpactRemoverWrapper{M<:MLJBase.Model} <: DeterministicComposite
	grp::Symbol
	classifier::M
	protected_attribute::Int64
	select_grp::Int64
end

function DisparateImpactRemoverWrapper(; classifier::MLJBase.Model=nothing, grp::Symbol=:class,protect_attribute::Int64, select_grp::Int64)
    model = DisparateImpactRemoverWrapper(grp, classifier, protected_attribute, select_grp)
    message = MLJBase.clean!(model)
    isempty(message)||@warn message
    return model
end

function MLJBase.clean!(model::DisparateImpactRemoverWrapper)
    warning = ""
    model.classifier!=nothing || (warning *= "No classifier specified in model\n")
    target_scitype(model) <: AbstractVector{<:Finite} || (warning *= "Only Binary Classifiers are supported\n")
    supports_weights(model.classifier) || (warning *= "Classifier provided does not support weights\n")
    return warning
end

function MLJBase.fit(model::DisparateImpactRemoverWrapper, verbosity::Int, X,y)
	grps = X[:, model.grp]

	X_repaired = DisparateImpactRemover(X, protected_attribute, select_grp)

	classifier = model.classifier
    mach1 = machine(classifier, X_repaired, y)
	fit!(mach1)
	return mach1
end

MMI.input_scitype(::Type{<:DisparateImpactRemoverWrapper{M}}) where M = input_scitype(M)
MMI.target_scitype(::Type{<:DisparateImpactRemoverWrapper{M}}) where M = AbstractVector{<:Finite{2}}
istype(model::DisparateImpactRemoverWrapper, type) = istype(model.classifier, type)


text:
using PyPlot: plot
## pyplot()

println("Ah) Generate and normalize continuum orbitals in a local potential.")

grid = Radial.Grid(Radial.Grid(false), rnt = 4.0e-6, h = 5.0e-2, hp = 2.0e-2, rbox = 20.0)
nuclearZ = 2.0;   noPoints = grid.NoPoints - 100

if  false
    # Define a full nuclear + Dirac-Fock-Slater potential
    wa   = Atomic.Computation(Atomic.Computation(), name="xx", grid=grid, nuclearModel=Nuclear.Model(nuclearZ),
                              configs=[Configuration("1s^2 2s"), Configuration("1s^2 2p")],
                              asfSettings=AsfSettings(AsfSettings(), scField=Basics.DFSField())  )

    wb = perform(wa, output=true)

    # Extract the ground level to generate a proper potential
    multiplet = wb["multiplet:"]
    ## grid      = wb["grid:"]
    level     = multiplet.levels[1]
    basis     = multiplet.levels[1].basis

    #
    println("\n\n>> Compute the nuclear and Hartree potentials.")
    w1  = JAC.Nuclear.pointNucleus(nuclearZ, grid)
    w2  = compute("radial potential: Dirac-Fock-Slater", grid, level)     ## compute("radial potential: Hartree",grid, level)
    ww = JAC.add(w1,w2)                                                   ## ww = Radial.Potential(" ", 0.95*ww.Zr, w2.grid)
else
    # Define a Coulomb potential
    ww  = JAC.Nuclear.pointNucleus(nuclearZ, grid)                         ## compute("radial potential: Hartree",grid, level)
    ## plot(grid.r, ww.Zr)
end
## plot(grid.r, ww.Zr);   error("stop aa")

if  true
    phases1 = Float64[];   phases2 = Float64[];   phases3 = Float64[];   phases4 = Float64[];
    norms1  = Float64[];   norms2  = Float64[];   norms3  = Float64[];   norms4  = Float64[]
    # Compute and normalize continuum orbitals for a range of energies
    settings = Continuum.Settings(false, noPoints)
    sh       = Subshell("1000s_1/2")
    energies = [x  for  x = 1.001:2.0:50.0];
    #
    @time for  energy in energies
        println("\n*** energy = $energy")
        wavenb   = sqrt( 2energy + energy * Defaults.getDefaults("alpha")^2 );      wavelgth = 2pi / wavenb
        println(">> Generate continuum orbital for kappa = $(sh.kappa) ($sh), energy = $energy a.u. and asympt. wavelength = $wavelgth; hp = $(grid.hp)")
        ## setDefaults("method: continuum, spherical Bessel");    setDefaults("method: normalization, pure sine")
        ## wc1      = Continuum.generateOrbitalLocalPotential(energy, sh, ww, settings)
        ## push!(phases1, wc1[2]);   push!(norms1, abs(wc1[3]))
        #
        setDefaults("method: continuum, Galerkin");            setDefaults("method: normalization, pure sine")
        wc2      = Continuum.generateOrbitalLocalPotential(energy, sh, ww, settings)
        push!(phases2, wc2[2]);   push!(norms2, abs(wc2[3]))
        setDefaults("method: continuum, Galerkin");            setDefaults("method: normalization, pure Coulomb")
        wc3      = Continuum.generateOrbitalLocalPotential(energy, sh, ww, settings)
        push!(phases3, wc3[2]);   push!(norms3, abs(wc3[3]))
        #
        setDefaults("method: continuum, Galerkin");            setDefaults("method: normalization, Ong-Russek")
        wc4      = Continuum.generateOrbitalLocalPotential(energy, sh, ww, settings)
        push!(phases4, wc4[2]);   push!(norms4, abs(wc4[3]))
    end
    #
    ##x plot(energies, norms2);   plot(energies, 1.005*norms3);   plot(energies, norms4) ## ;   plot(energies, phases4)
    ##x plot(energies, phases2);   plot(energies, phases3);   plot(energies, phases4);   plot(energies, norms2)
    plot(energies, phases2);   plot(energies, 1.005*phases3);   plot(energies, 1.01*phases4)

else
    # Compute and display single continuum orbital
    ##  "method: continuum, pure sine"  "method: continuum, asymptotic Coulomb"  "method: continuum, Galerkin"  "method: continuum, spherical Bessel"
    ##  "method: normalization, pure sine"    "method: normalization, pure Coulomb"  "method: normalization, Ong-Russek"
    @show Defaults.GBL_CONT_SOLUTION, Defaults.GBL_CONT_NORMALIZATION
    settings = Continuum.Settings(false, noPoints)
    sh       = Subshell("1000g_7/2")
    energies = [0.5];
    for  energy in energies
        wavenb   = sqrt( 2energy + energy * Defaults.getDefaults("alpha")^2 );      wavelgth = 2pi / wavenb
        println(">> Generate continuum orbital for kappa = $(sh.kappa) ($sh), energy = $energy a.u. and asympt. wavelength = $wavelgth; hp = $(grid.hp)")
        @show " "
        ## setDefaults("method: continuum, pure sine");           setDefaults("method: normalization, pure sine")
        ## @show Defaults.GBL_CONT_SOLUTION, Defaults.GBL_CONT_NORMALIZATION
        ## wc1      = Continuum.generateOrbitalLocalPotential(energy, sh, ww, settings)
        ## plot(grid.r[1:noPoints], wc1[1].P[1:noPoints])
        #
        ## setDefaults("method: continuum, spherical Bessel");    setDefaults("method: normalization, pure sine")
        ## @show " "
        ## @show Defaults.GBL_CONT_SOLUTION, Defaults.GBL_CONT_NORMALIZATION
        ## wc2      = Continuum.generateOrbitalLocalPotential(energy, sh, ww, settings)
        ## plot(grid.r[1:noPoints], wc2[1].P[1:noPoints])
        #
        setDefaults("method: continuum, Galerkin");            setDefaults("method: normalization, pure sine")
        @show " "
        @show Defaults.GBL_CONT_SOLUTION, Defaults.GBL_CONT_NORMALIZATION
        wc3      = Continuum.generateOrbitalLocalPotential(energy, sh, ww, settings)
        plot(grid.r[1:noPoints], wc3[1].P[1:noPoints])
        #
        setDefaults("method: continuum, Galerkin");            setDefaults("method: normalization, pure Coulomb")
        @show " "
        @show Defaults.GBL_CONT_SOLUTION, Defaults.GBL_CONT_NORMALIZATION
        wc4      = Continuum.generateOrbitalLocalPotential(energy, sh, ww, settings)
        plot(grid.r[1:noPoints], wc4[1].P[1:noPoints])
        #
        setDefaults("method: continuum, Galerkin");            setDefaults("method: normalization, Ong-Russek")
        @show " "
        @show Defaults.GBL_CONT_SOLUTION, Defaults.GBL_CONT_NORMALIZATION
        wc5      = Continuum.generateOrbitalLocalPotential(energy, sh, ww, settings)
        plot(grid.r[1:noPoints], wc5[1].P[1:noPoints])
    end
end


text:
################################################################################
#
#  Algorithm of Schmettow
#
################################################################################

function one_step(b::NfOrdFracIdl, p::NfOrdIdl; prec::Int = 100)
  b = p*b
  simplify(b)
  g1 = short_elem(b, prec = prec)
  b = g1*inv(b)
  simplify(b)
  g2 = short_elem(b, prec = prec)
  return simplify(g2*inv(b)), g1, g2
end


text:
module TypingTests

using Interpolations, Base.Test

nx = 10
f(x) = convert(Float32, x^3/(nx-1))
g(x) = convert(Float32, 3x^2/(nx-1))

A = Float32[f(x) for x in 1:nx]

itp = interpolate(A, BSpline(Quadratic(Flat())), OnCell())

# display(plot(
#     layer(x=1:nx,y=[f(x) for x in 1:1//1:nx],Geom.point),
#     layer(x=1:.1:nx,y=[itp[x] for x in 1:1//10:nx],Geom.path),
# ))

for x in 3.1:.2:4.3
    @test_approx_eq_eps float(f(x)) float(itp[x]) abs(.1*f(x))
end

@test typeof(itp[3.5f0]) == Float32

for x in 3.1:.2:4.3
    @test_approx_eq_eps g(x) gradient(itp, x) abs(.1*g(x))
end

@test typeof(gradient(itp, 3.5f0)[1]) == Float32

# Rational element types
R = Rational{Int}[x^2//10 for x in 1:10]
itp = interpolate(R, BSpline(Quadratic(Free())), OnCell())
itp[11//10]

@test typeof(itp[11//10]) == Rational{Int}
@test itp[11//10] == (11//10)^2//10

end


text: circularPermutation(n, start) = map(x -> start ⊻ x ⊻ (x >> 1) ,0:((1<<n) - 1) )

max_stars_repo_name: ccqpein/Arithmetic-Exercises

max_issues_repo_name: ccqpein/Arithmetic-Exercises

max_forks_repo_name: ccqpein/Arithmetic-Exercises


text:
"""
    Molmer_QFI_GHZ(Nj, Tfinal, dt; kwargs...)

Numerically exact evaluation of the Mølmer QFI for the GHZ state for ``Nj`` atoms interacting with independent noise. `Tfinal` and `dt` specify the

Returns `(t, qfi)` where `qfi` is the Mølmer QFI evaluated at time instants in `t`.

# Keyword arguments

    * `κ = 1`: coupling with the noise
    * `θ = 0`: angle with the noise
    * `ω = 0`: magnetic field / frequency
    * `dω = 1e-5`: increment dω for numerically evaluating the derivative
"""
function Molmer_QFI_GHZ(
    Nj::Int64,      # Number of spins
    Tfinal,         # Final time
    dt::Float64;    # Time step
    κ = 1.,         # Field rate
    θ = 0.,         # Angle of the noise
    ω = 0.,         # Magnetic field
    dω=1e-5        # Increment to numerically compute the derivative (instable if too small wrt to ω)
    )

    Ntime = Int(floor(Tfinal/dt))
    t = (1 : Ntime) * dt
    θ = pi/2 - θ # Matteo chose the opposite convention for the effective QFI
    QFisherT = zero(t)

    lmat=[ 0 0 0 -1im.*dω ;
            0  -κ.*sin(θ).^2  -ω   κ.*cos(θ).*sin(θ) ;
            0 ω -κ 0 ;
             -1im.*dω  κ.*cos(θ).*sin(θ) 0 -κ.*cos(θ).^2]

    for jt=1:Ntime
        fmat=expm(lmat*t[jt])
        f00=sqrt(0.5).*fmat*[1. 0. 0. 1.]'
        f11=sqrt(0.5).*fmat*[1. 0. 0. -1.]'
        f01=sqrt(0.5).*fmat*[0. 1. 1im 0.]'
        f10=sqrt(0.5).*fmat*[0. 1. -1im 0.]'
        logtrρ=log.(abs(0.5 .*( (f00[1] .* sqrt(2.)) .^ Nj + (f11[1] .* sqrt(2.)) .^ Nj + (f01[1] .* sqrt(2.)) .^ Nj + (f10[1] .* sqrt(2.)) .^ Nj )))
        # check: the trace of ρ should be less then 1... it can be approximately 0, but positive:
        # in that case we manually set it to 0
        if logtrρ<0
            QFisherT[jt] = 2. .* exp.(log.(-logtrρ) - 2. *log.(dω))
        else
            QFisherT[jt] = 0.
        end
    end
    return (t,QFisherT)
end

"""
    Molmer_qfi_transverse(t, Nj, κ)

The Mølmer QFI for ``N_j`` atoms interacting with independent transverse
noise with coupling ``\\κ``.

`t` can be any `Array`.

This function effectively evaluates the formula

``Q(t) = 4  \\frac{N_j}{\\κ ^2} e^{-\\κ  t} \\left(2 (2- N_j) e^{\\frac{\\κ  t}{2}}+e^{\\κ
   t} (\\κ  t+ N_j-3)+N_j-1\\right)``
"""
function Molmer_qfi_transverse(t, Nj, κ)
    return exp.(-2 * κ .* t) .* κ.^(-2) .* Nj  .*
     (-1 -2 .* exp.(κ .* t) .* (Nj - 2) + Nj +
        exp.(2. * κ .* t ) .* (-3+Nj +2. * κ .* t));
end


text:
const PROCESSING_STEPS = (
    Predict = LinearModel,
    Cluster = Cluster,
    Project = DimensionalityReduction,
    Wildcard = Wildcard,
)

struct Process <: AbstractPipeline
    table::Observable{SimpleTable}
    list::EditableList
    value::Observable{SimpleTable}
    options::Vector{Symbol}
end

get_vertices(p::Process) = Vertex.(p.list.steps[])
get_vertex_names(p::Process) = p.options

function Process(table::Observable{SimpleTable}, options::Vector{Symbol})
    value = Observable(table[])
    steps = Observable(Any[])
    function thunkify(type)
        return function ()
            step = type(value)
            card = step.card
            on(card.run) do _
                if shouldrun(card.state[])
                    value[] = compute_pipeline(table[], value[], steps[])
                end
            end
            on(card.destroy) do val
                clear!(card)
                _steps = steps[]
                if val
                    id = objectid(step)
                    idx = findfirst(==(id)∘objectid, _steps)
                    isnothing(idx) || (steps[] = remove_item(_steps, idx))
                end
            end
            return step
        end
    end
    _list_options = Any[SimpleDict("key" => "Move Selected", "value" => MoveSelected)]
    for option in options
        option_dict = SimpleDict(
            "key" => string(option),
            "value" => thunkify(PROCESSING_STEPS[option])
        )
        push!(_list_options, option_dict)
    end
    list_options = Observable(_list_options)
    process = Process(table, EditableList(list_options, steps), value, options)
    on(table) do data
        value[] = compute_pipeline(always_true, data, value[], steps[])
    end
    return process
end

function Process(table::Observable{SimpleTable}; options=[:Predict, :Cluster, :Project])
    return Process(table, vecmap(Symbol, options))
end

function jsrender(session::Session, process::Process)
    ui = scrollable_component(
        process.list;
        onmousedown=js"$(UtilitiesJS).updateSelection(this, event, $(process.list.selected));"
    )
    return jsrender(session, ui)
end

default_needs_update(step) = shouldrun(step.card.state[])
always_true(step) = true

Vertex(step::AbstractProcessingStep) = Vertex(step.card.name, columns_in(step), columns_out(step))

nodes_to_compute(steps) = nodes_to_compute(default_needs_update, steps)

function nodes_to_compute(f, steps)
    g = simpledigraph(Vertex.(steps))
    sorted = topological_sort_by_dfs(g)
    needs_updating = map(f, steps)
    nodes = Int[]
    for node in sorted
        for neighbor in inneighbors(g, node)
            needs_updating[node] |= needs_updating[neighbor]
        end
        needs_updating[node] && push!(nodes, node)
    end
    return nodes
end

compute_pipeline(input, cache, steps) = compute_pipeline(default_needs_update, input, cache, steps)

function compute_pipeline(f, input, cache, steps)
    nodes =  nodes_to_compute(f, steps)
    result = SimpleTable(input)
    # Add outputs of nodes that do not need to be recomputed
    for node in setdiff(1:length(steps), nodes)
        names = columns_out(steps[node])
        columns = AbstractVector[Tables.getcolumn(cache, name) for name in names]
        mergedisjointcols!(result, SimpleTable(names, columns))
    end
    # Recompute nodes, which were returned in topological order
    for node in nodes
        step = steps[node]
        if !isempty(columns_in(step))
            try
                mergedisjointcols!(result, step(result))
                step.card.state[] = done
            catch e
                step.card.error[] = sprint(showerror, e)
                step.card.state[] = errored
            end
        else
            step.card.state[] = inactive
        end
    end
    return result
end


text:
using ResponseManneville
using Test

@testset "ResponseManneville.jl" begin

include("testLSV.jl")
include("testC2Basis.jl")

end

max_stars_repo_name: orkolorko/ResponseManneville

max_issues_repo_name: orkolorko/ResponseManneville

max_forks_repo_name: orkolorko/ResponseManneville


text:
const typecodes = Dict{Int32, DataType}([
    0 => UInt8,
    1 => Float32,
    2 => Float64,
    3 => Int32
])

using SparseArrays

function deserialize(fname::AbstractString)
    file = open(fname, "r")
    typecode = read(file, Int32)
    if !(typecode in keys(typecodes))
        throw(ErrorException("Unknown typecode"))
    end

    dtype = typecodes[typecode]

    I = UInt64[]
    J = UInt64[]
    V = dtype[]

    nrows = read(file, UInt64)
    row_lens = resize!(UInt64[], nrows)
    read!(file, row_lens)

    nentries = zero(UInt64)
    for i = 1:nrows
        nentries += row_lens[i]
        append!(I, (i for j = 1:row_lens[i]))
    end

    resize!(J, nentries)
    resize!(V, nentries)
    read!(file, J)
    read!(file, V)

    sparse(I, J .+= 1, V)
end

max_stars_repo_name: slmcbane/SymSparse

max_issues_repo_name: slmcbane/SymSparse

max_forks_repo_name: slmcbane/SymSparse


text:
# Keep this code as a library for future use. The JacobiSVD is significantly faster than the SVD that comes standard with Julia.

function ATA(A::AbstractArray{T,2}, j::Int, k::Int) where {T<:Real} # need AbstractArray instead of Array to handle Transpose, which is an AbstractArray but not an Array
    sum = zero(T)
    nrows, _ = size(A)
    @simd for i in 1:nrows
        sum += A[i, j] * A[i, k]
    end
    return sum
end

@inline epsilon(T::Real) = 4 * eps(T)

"this version will always perform iterations to high relative accuracy"
GivensRotation(ajj::T, ajk::T, akk::T) where {T<:Real} = GivensRotation(false, ajj, ajk, akk)

function GivensRotation(lowRelativeAccuracy::Bool, ajj::T, ajk::T, akk::T) where {T<:Real}
    local sinTheta, cosTheta, needToRotate

    product = ajj * akk # default value that may be overriden if computing small singular values to low relative accuracy

    if lowRelativeAccuracy
        println("using low relative accuracy")
        # less stringent accuracy requirement than requiring that the rotation matrices not change the diagonal elements of the matrix at all. The result is that small singular values may be computed with poor
        # relative accuracy, possibly having few significant digits of accuracy. For most applications computing the small singular values to high accuracy is not important so this is probably the way you would
        # want to compute most of the time.
        max = max(abs(ajj), abs(akk))
        min = min(abs(ajj), abs(akk))
        if min == 0 || max / min > 2^52  # this will cause the test "if (abs(ajk) >= eps * sqrt(product))" to fail with fewer iterations which makes the code significantly more efficient. This needs to be rewritten in terms of machine epsilon. This number pow(2,52) is correct only for Float64.
            product = max
        end
    end

    if ajk != zero(T) && abs(ajk) >= eps(T) * sqrt(product)
        local c, s, t, tau

        tau = (ajj - akk) / (2 * ajk)
        t = sign(tau) / (abs(tau) + sqrt(1 + tau * tau))
        c = 1 / sqrt(1 + t * t)
        s = c * t
        sinTheta = s
        cosTheta = c

        # NOTE: commented out this code because it greatly increased the residual Ax-b. Not sure how to include this without messing up accuracy in general. Leave it out for now.

        # when the matrix is singular it can be the case that the diagonal elements ajj, akk go to zero faster than ajk. This means the test "if (abs(ajk) >= eps * sqrt(ajj*akk))"
        # returns true until the diagonal and off diagonal elements become so small that the rotation matrix has no effect on any of the elements in the matrix (the rotation matrix becomes almost exactly
        # equal to the identity). If this happens the JacobiSVD function
        # will never stop iterating and you will get a max iterations exceeded exception. The following test detects when the diagonal elements are no longer being changed by rotation
        # and causes the iteration to stop if this is the case.
        # if abs((c * ajj - s * ajk) - ajj) < eps(T) && abs((s * ajk + c * akk) - akk) < eps(T)
        #     # c and s values are such that none of the diagonal elements in the matrix will change if the rotation is performed so don't rotate. More iterations than necessary since have to wait till diagonal and
        #     # off diagonal elements are reduced to values so small they are beyond the range of whatever precision float you are using. This is probably necessary if want high relative accuracy in small singular values but
        #     # there are many cases when this is not a requirement.
        #     # Have to use eps(T) instead of a test for equality because C# is apparently not using 64 bit IEEE double precision arithmetic. Instead it appears to be using the Intel extended precision
        #     # floating point. It takes many more iterations to reduce the diagonal elements to zero in this precision than in double precision.
        #     # an even more stringent requirement would be for the off diagonal elements also not to be changed by the rotation matrix but this can increase the number of iterations by 3-5 times for small
        #     # matrices, maybe more for big matrices. In some extreme cases exiting before the off diagonal elements cease to be changed may result in more error in the singular vectors than desirable but I have not observed this yet.
        #     needToRotate = false
        #     return sinTheta,cosTheta,needToRotate
        # else
        #     needToRotate = true
        #     return sinTheta,cosTheta,needToRotate
        # end

        needToRotate = true
        return sinTheta, cosTheta, needToRotate
    else
        # off diagonal elements essentially zero
        sinTheta = 0
        cosTheta = 0
        needToRotate = false
        return sinTheta, cosTheta, needToRotate
    end
end

function oneSidedJacobiRotation!(A::AbstractArray{T,2}, j::Int, k::Int, sinTheta::T, cosTheta::T) where {T<:Real} # need AbstractArray instead of Array to handle Transpose, which is an AbstractArray but not an Array
    s = sinTheta
    c = cosTheta
    numrows, _ = size(A)

    for i in 1:numrows
        Aij = A[i, j]
        Aik = A[i, k]
        A[i, j] = c * Aij + s * Aik
        A[i, k] = -s * Aij + c * Aik
    end
end

"use this version if you want to handle exceeding maximum number of iterations. DESTRUCTIVE: matrix A is overwritten. powerPercentage must be be between 0 and 1"
function JacobiSolve(A::AbstractArray{T,2}, b::AbstractArray{T,1}; computeToLowRelativeAccuracy::Bool = false, conditionNumber::T = T(1e20))::Array{T,1} where {T<:Real}  # don't like having to declare this as Real rather than AbstractFloat because eps(Real) wouldn't necessarily be defined. But that's the way the IntervalArithmetic package is designed so to make it work with intervals have to this. The ForwardDiff package is also designed to take Real arguments.
    local temp, total
    currentSum = zero(T)

    U, sigma, V, status = JacobiSVD(computeToLowRelativeAccuracy, A)
    temp = U' * b

    for i in 1:length(sigma)
        # println(S[0]/S[S.GetLength(0)-1]);
        if sigma[1] / sigma[i] > conditionNumber || sigma[i] == zero(T)
            temp[i] = zero(T)
        else
            temp[i] /= sigma[i]
        end
    end

    return V * temp

end

"DESTRUCTIVE: matrix A is overwritten. Returns U,sigma,V (not V'), iterationsexceeded   where iterationsexceeded will be either :exceededMax or :didNotExceedMax. If max iterations was exceeded it is possible the result is not accurate. But the results might also be fine."
function JacobiSVD(relativeaccuracy::Bool, A::AbstractArray{T,2}) where {T<:Real}
    acopy = copy(A)
    return JacobiSVD!(relativeaccuracy, acopy)
end

"sets flag to compute singular values to high relative accuracy. DESTRUCTIVE: matrix A is overwritten"
function JacobiSVD!(A::AbstractArray{T,2}) where {T<:Real}
    # out double[,] U, out double[] sigma, out double[,] V)
    return JacobiSVD!(false, A)
end

"DESTRUCTIVE: matrix A is overwritten. Returns U,sigma,V (not V'), iterationsexceeded   where iterationsexceeded will be either :exceededMax or :didNotExceedMax. If max iterations was exceeded it is possible the result is not accurate. But the results might also be fine."
function JacobiSVD!(computeToLowRelativeAccuracy::Bool, A::AbstractArray{T,2}) where {T<:Real}
    local status

    maxIterations = 35
    computingTranspose = false
    rowsA, colsA = size(A)

    # extend later
    if rowsA < colsA
        A = transpose(A)
        computingTranspose = true
        rowsA, colsA = size(A) # compute new size of transposed A
    end

    nrows = rowsA
    ncols = colsA
    n = colsA
    minRowsCols = min(nrows, ncols)
    sigma = Array{T,1}(undef, minRowsCols)
    localSigma = similar(sigma)

    numIterations = 0
    R = Matrix{T}(I, ncols, ncols)

    V = Array{T,2}(undef, n, n)
    U = Array{T,2}(undef, nrows, minRowsCols)
    while true

        rotatedThisIteration = false

        for j in 1:(n - 1)
            for k in (j + 1):n
                GTGjj = ATA(A, j, j)
                GTGjk = ATA(A, j, k)
                GTGkk = ATA(A, k, k)
                s, c, needToRotate = GivensRotation(computeToLowRelativeAccuracy, GTGjj, GTGjk, GTGkk)
                rotatedThisIteration = rotatedThisIteration || needToRotate
                if needToRotate
                    oneSidedJacobiRotation!(A, j, k, s, c)
                    oneSidedJacobiRotation!(R, j, k, s, c)
                end
            end
        end

        numIterations += 1
        if numIterations > maxIterations
            break
        end
        if !rotatedThisIteration
            break
        end
    end

    for i in 1:length(localSigma)
        localSigma[i] = 0
        for j in 1:rowsA
            # localSigma[i] += A[j, i] * A[j, i];
            localSigma[i] += A[j, i]^2 # use this instead of A[j,i]*A[j,i]. Interval arithmetic package will compute a much tighter inclusion
        end
        localSigma[i] = sqrt(localSigma[i])
    end

    # sort sigmas
    indices = sortperm(localSigma, rev = true)
    localSigma = localSigma[indices]

    for i in 1:minRowsCols
        sigma[i] = localSigma[i]  # copy just the true singular values

        for j in 1:rowsA
            if (sigma[i] != 0)
                U[j, i] = A[j, indices[i]] / sigma[i]
            else
                U[j, i] = 0
            end

        end
        for j in 1:n
            V[j, i] = R[j, indices[i]]
        end
    end

    if computingTranspose
        temp = V
        V = U
        U = temp
    end

    if numIterations > maxIterations
        status = :exceededMax
    else
        status = :didNotExceedMax
    end
    return (U = U, sigma = sigma, V = V, status = status)
end


text:
module JACKAudioTests

using Compat
import Compat.ASCIIString

if VERSION >= v"0.5.0-"
    using Base.Test
else
    using BaseTestNext
end
using JACKAudio
using SampledSignals

function jackver()
    verlines = readlines(ignorestatus(`jackd --version`))
    println(join(verlines))
    jackd = split(verlines[1])[1]

    jackd == "jackdmp" ? 2 : 1
end

const jackd = jackver()

# GENERAL WARNING: you probably want to run jackd with as large a buffer as
# possible (e.g. 2048) to reduce the chance of nondeterminism caused by the
# jack callback happening when you don't want it to

@testset "JACK Tests" begin
    # we call the read and write functions here to make sure they're precompiled
    c = JACKClient()
    buf = read(c, 1)
    read!(c, buf)
    write(c, buf)
    close(c)

    c = JACKClient(1, 1)
    buf = read(c, 1)
    read!(c, buf)
    write(c, buf)
    # and a 1D buf
    monobuf = SampleBuf(rand(Float32, 32), samplerate(c))
    read!(c, monobuf)
    write(c, monobuf)
    close(c)

    # the process callback is not part of the public API, but we want to run
    # some tests on it anyways. This seems to segfault on jack1
    println("test 1")
    jackd == 2 && @testset "Process Callback" begin
        client = JACKClient(active=false)
        # note we're caching the client.portptrs access because it seems to
        # cause 16 bytes of allocation
        ptrs = client.portptrs
        # make sure we run it to warm up
        JACKAudio.process(UInt32(256), ptrs)
        alloc = @allocated JACKAudio.process(UInt32(256), ptrs)
        close(client)
        @test alloc == 0
    end
    println("test 2")
    @testset "No-Argument Construction" begin
        c = JACKClient()
        @test c.name == "Julia"
        @test length(sources(c)) == 1
        @test nchannels(sources(c)[1]) == 2
        @test length(sinks(c)) == 1
        @test nchannels(sinks(c)[1]) == 2
        close(c)
    end
    println("test 3")
    @testset "Channel Count Construction" begin
        c = JACKClient(4, 5)
        @test c.name == "Julia"
        @test length(sources(c)) == 1
        @test nchannels(sources(c)[1]) == 4
        @test length(sinks(c)) == 1
        @test nchannels(sinks(c)[1]) == 5
        close(c)
    end
    println("test 4")
    @testset "Name Construction" begin
        c = JACKClient("TestClient")
        @test c.name == "TestClient"
        @test length(sources(c)) == 1
        @test nchannels(sources(c)[1]) == 2
        @test length(sinks(c)) == 1
        @test nchannels(sinks(c)[1]) == 2
        close(c)
    end
    println("test 5")
    @testset "Full Custom Construction" begin
        c = JACKClient("TestClient", [("In1", 2), ("In2", 3)],
                                     [("Out1", 1), ("Out2", 2)])
        @test c.name == "TestClient"
        @test length(sources(c)) == 2
        @test nchannels(sources(c)[1]) == 2
        @test nchannels(sources(c)[2]) == 3
        @test length(sinks(c)) == 2
        @test nchannels(sinks(c)[1]) == 1
        @test nchannels(sinks(c)[2]) == 2
        close(c)
    end
    println("test 6")
    # this test is an example of how finicky it is to do synchronized audio IO
    # using a stream-based read/write API.
    @testset "Short Reading/Writing (less than ringbuffer) works" begin
        sourceclient = JACKClient("Source", 1, 0; connect=false)
        sinkclient = JACKClient("Sink", 0, 1; connect=false)
        # connect them in JACK
        connect(sinkclient, sourceclient)
        buf = SampleBuf(rand(Float32, 32, 1), samplerate(sourceclient))
        readbuf = SampleBuf(Float32, samplerate(sourceclient), size(buf))
        # clear out any audio we've accumulated in the source buffer
        seekavailable(sources(sourceclient)[1])
        # synchronize so we know we're running the test at the beginning of the
        # block. It's still not 100% deterministic but hopefully this makes it
        # pass most of the time.
        read(sourceclient, 1)
        # skip the rest of the block we received
        seekavailable(sources(sourceclient)[1])
        # this won't block because we just write directly into the ringbuf
        write(sinkclient, buf)
        # this should block now as well because there weren't any more frames
        # to read. In the next process callback JACK should read from the sink
        # and write to the source, sticking the data in readbuf
        read!(sourceclient, readbuf)
        @test buf == readbuf
        close(sourceclient)
        close(sinkclient)
    end
    println("test 7")

    @testset "Long reading/writing (more than ringbuffer) works" begin
        sourceclient = JACKClient("Source", 1, 0; connect=false)
        sinkclient = JACKClient("Sink", 0, 1; connect=false)
        connect(sinkclient, sourceclient)

        writebuf = SampleBuf(rand(Float32, JACKAudio.RINGBUF_SAMPLES + 32), samplerate(sourceclient))
        readbuf = SampleBuf(Float32, samplerate(sourceclient), size(writebuf))

        seekavailable(sources(sourceclient)[1])
        read(sourceclient, 1)
        seekavailable(sources(sourceclient)[1])
        @sync begin
            @async write(sinkclient, writebuf)
            @async read!(sourceclient, readbuf)
        end
        @test readbuf == writebuf
        close(sourceclient)
        close(sinkclient)
    end
    println("test 8")

    @testset "writers get queued" begin
        sourceclient = JACKClient("Source", 1, 0; connect=false)
        sinkclient = JACKClient("Sink", 0, 1; connect=false)
        connect(sinkclient, sourceclient)

        writebuf1 = SampleBuf(rand(Float32, JACKAudio.RINGBUF_SAMPLES * 2 + 32), samplerate(sourceclient))
        writebuf2 = SampleBuf(rand(Float32, JACKAudio.RINGBUF_SAMPLES * 2 + 32), samplerate(sourceclient))
        readbuf = SampleBuf(Float32, samplerate(sourceclient), length(writebuf1) + length(writebuf2))

        seekavailable(sources(sourceclient)[1])
        read(sourceclient, 1)
        seekavailable(sources(sourceclient)[1])
        @sync begin
            @async write(sinkclient, writebuf1)
            @async write(sinkclient, writebuf2)
            @async read!(sourceclient, readbuf)
        end
        @test readbuf[1:length(writebuf1)] == writebuf1
        @test readbuf[(length(writebuf1)+1):end] == writebuf2
        close(sourceclient)
        close(sinkclient)
    end

    println("test 9")

    @testset "readers get queued" begin
        sourceclient = JACKClient("Source", 1, 0; connect=false)
        sinkclient = JACKClient("Sink", 0, 1; connect=false)
        connect(sinkclient, sourceclient)

        writebuf1 = SampleBuf(rand(Float32, JACKAudio.RINGBUF_SAMPLES * 2 + 32), samplerate(sourceclient))
        writebuf2 = SampleBuf(rand(Float32, JACKAudio.RINGBUF_SAMPLES * 2 + 32), samplerate(sourceclient))
        readbuf = SampleBuf(Float32, samplerate(sourceclient), length(writebuf1) + length(writebuf2))

        seekavailable(sources(sourceclient)[1])
        read(sourceclient, 1)
        seekavailable(sources(sourceclient)[1])
        @sync begin
            @async write(sinkclient, writebuf1)
            @async write(sinkclient, writebuf2)
            @async read!(sourceclient, readbuf)
        end
        @test readbuf[1:length(writebuf1)] == writebuf1
        @test readbuf[(length(writebuf1)+1):end] == writebuf2
        close(sourceclient)
        close(sinkclient)
    end
end

end # module