text:
"""
    mutable struct AWSConfig
        account_id::Union{Missing, String} = missing
        region::Union{Missing, String} = missing
    end

Defines the Amazon Web Services account id and region to use.

Your current AWS CLI profile should match the account id specified here. To check that this is the
case, run `aws sts get-caller-identity --query Account --output text` from the command line. To
see available profiles, run `aws configure list-profiles`.
"""
@with_kw mutable struct AWSConfig
  account_id::Union{Missing, String} = missing
  region::Union{Missing, String} = missing
end
StructTypes.StructType(::Type{AWSConfig}) = StructTypes.Mutable()

function get_aws_config()::AWSConfig
  AWSConfig(
            readchomp(`aws sts get-caller-identity --query Account --output text`),
            readchomp(`aws configure get region`)
           )
end

@with_kw mutable struct AWSRolePolicyStatement
  Effect::Union{Missing, String} = missing
  Principal::Union{Missing, Dict{String, Any}} = missing
  Action::Union{Missing, String} = missing
end
StructTypes.StructType(::Type{AWSRolePolicyStatement}) = StructTypes.Mutable()
Base.:(==)(a::AWSRolePolicyStatement, b::AWSRolePolicyStatement) = (
  a.Effect == b.Effect && a.Principal == b.Principal && a.Action == b.Action)

@with_kw mutable struct AWSRolePolicyDocument
  Version::Union{Missing, String} = missing
  Statement::Vector{AWSRolePolicyStatement} = Vector{AWSRolePolicyStatement}()
end
StructTypes.StructType(::Type{AWSRolePolicyDocument}) = StructTypes.Mutable()
Base.:(==)(a::AWSRolePolicyDocument, b::AWSRolePolicyDocument) = (a.Version == b.Version && a.Statement == b.Statement)

@with_kw mutable struct AWSRole
  Path::Union{Missing, String} = missing
  RoleName::Union{Missing, String} = missing
  RoleId::Union{Missing, String} = missing
  Arn::Union{Missing, String} = missing
  CreateDate::Union{Missing, String} = missing
  AssumeRolePolicyDocument::Union{Missing, AWSRolePolicyDocument} = missing
  MaxSessionDuration::Union{Missing, Int64} = missing
  exists::Bool = true
end
StructTypes.StructType(::Type{AWSRole}) = StructTypes.Mutable()
Base.:(==)(a::AWSRole, b::AWSRole) = a.RoleId == b.RoleId

const lambda_execution_policy_statement = AWSRolePolicyStatement(
    Effect = "Allow",
    Principal = Dict("Service" => "lambda.amazonaws.com"),
    Action = "sts:AssumeRole",
  )

"""
    function get_all_aws_roles()::Vector{AWSRole}
Get all AWS roles from the current AWS Config.
"""
function get_all_aws_roles()::Vector{AWSRole}
  all_roles_json = readchomp(`aws iam list-roles`)
  all = JSON3.read(all_roles_json, Dict{String, Vector{AWSRole}})
  all["Roles"]
end

"""
    function get_aws_role(role_name::String)::Union{Nothing, AWSRole}
Get the AWS Role corresponding to the given role_name. Returns nothing if none found.
"""
function get_aws_role(role_name::String)::Union{Nothing, AWSRole}
  all = get_all_aws_roles()
  index = findfirst(role -> role.RoleName == role_name, all)
  isnothing(index) ? nothing : all[index]
end

"""
    function create_aws_role(role_name::String)::AWSRole

Create an AWS Role with Lambda execution permissions.
"""
function create_aws_role(role_name::String)::AWSRole
  existing_role = get_aws_role(role_name)
  @debug role_name
  @debug existing_role
  if !isnothing(existing_role)
    if aws_role_has_lambda_execution_permissions(existing_role)
      existing_role
    else
      error("Existing role with name $role_name exists, but does not have execution permissions")
    end
  else
    create_script = get_create_lambda_role_script(role_name)
    role_json = readchomp(`bash -c $create_script`)
    @info "Creating role $role_name ..."
    sleep(8);
    JSON3.read(role_json, Dict{String, AWSRole})["Role"]
  end
end

"""
    function delete!(role::AWSRole)

Delete an AWS Role. The role itself will be deleted from AWS, and the `AWSRole` type will no longer
be usable.
"""
function delete!(role::AWSRole)
  role.exists || error("Role does not exist")
  delete_script = get_delete_lambda_role_script(role.RoleName)
  run(`bash -c $delete_script`)
  role.exists = false
end

function get_role_arn_string(
    aws_config::AWSConfig,
    role_name::String,
  )::String
  "arn:aws:iam::$(aws_config.account_id):role/$role_name"
end

function get_function_uri_string(aws_config::AWSConfig, function_name::String)::String
  "$(aws_config.account_id).dkr.ecr.$(aws_config.region).amazonaws.com/$function_name"
end

function get_function_arn_string(aws_config::AWSConfig, function_name::String)::String
  "arn:aws:lambda:$(aws_config.region):$(aws_config.account_id):function:$function_name"
end

function get_ecr_arn_string(aws_config::AWSConfig, image_suffix::String)::String
  "arn:aws:ecr:$(aws_config.region):$(aws_config.account_id):repository/$image_suffix"
end

function get_ecr_uri_string(aws_config::AWSConfig, image_suffix::String)::String
  "$(aws_config.account_id).dkr.ecr.$(aws_config.region).amazonaws.com/$image_suffix"
end

function create_lambda_execution_role(role_name)
  create_script = get_create_lambda_role_script(role_name)
  run(`bash -c $create_script`)
end

function aws_role_has_lambda_execution_permissions(role::AWSRole)::Bool
  lambda_execution_policy_statement in role.AssumeRolePolicyDocument.Statement
end


text:
using Test
using LinearAlgebra
using PolynomialBasis
using ImplicitDomainQuadrature
# using Revise
using HDGElasticity

function allapprox(v1,v2)
    return all(v1 .≈ v2)
end

sbasis = TensorProductBasis(1,1)
squad = tensor_product_quadrature(1,2)
facequads = repeat([squad],4)

testmatrix = [2/3 0.0 1/3 0.0
              0.0 2/3 0.0 1/3
              1/3 0.0 2/3 0.0
              0.0 1/3 0.0 2/3]

HH = HDGElasticity.HHop(sbasis,squad,1.)
@test allapprox(HH,testmatrix)

cellmap = HDGElasticity.CellMap([0.,0.],[1.,0.5])
facescale = HDGElasticity.face_determinant_jacobian(cellmap)
HH = [HDGElasticity.HHop(sbasis,facequads[i],facescale[i]) for i = 1:4]

testmatrix = [2/3 0.0 1/3 0.0
              0.0 2/3 0.0 1/3
              1/3 0.0 2/3 0.0
              0.0 1/3 0.0 2/3]

@test allapprox(HH[1],0.5*testmatrix)
@test allapprox(HH[2],0.25*testmatrix)
@test allapprox(HH[3],0.5*testmatrix)
@test allapprox(HH[4],0.25*testmatrix)

HH = reshape([0.0],1,1)
sfunc(x) = [0.5*x[1]*(x[1]+1)]
squad = tensor_product_quadrature(1,3)
components = [1.]

HDGElasticity.HHop!(HH,sfunc,squad,components,1.,1)
@test allapprox(HH,[4/15])

fill!(HH,0.0)
HDGElasticity.HHop!(HH,sfunc,squad,components,2.,1)
@test allapprox(HH,2*[4/15])

fill!(HH,0.0)
components = [1.0]
HDGElasticity.HHop!(HH,sfunc,squad,components,3.,1)
@test allapprox(HH,3*[4/15])

n = [1.,1.]/sqrt(2.)
normals = repeat(n,inner=(1,3))
imap = InterpolatingPolynomial(2,sbasis)
coeffs = [-0.5  -1.0
          -1.0  +0.0]
update!(imap,coeffs)
cellmap = HDGElasticity.CellMap([0.,0.],[2.,1.])

HH = zeros(4,4)
facescale = HDGElasticity.scale_area(cellmap,normals)
HDGElasticity.HHop!(HH,sbasis,squad,imap,facescale,2)
@test allapprox(HH,0.5/sqrt(2.0)*testmatrix)

fill!(HH,0.0)
HDGElasticity.HHop!(HH,sbasis,squad,imap,2facescale,2)
@test allapprox(HH,1.0/sqrt(2.0)*testmatrix)

HH = HDGElasticity.HHop(sbasis,squad,imap,normals,cellmap)
@test allapprox(HH,0.5/sqrt(2.0)*testmatrix)

cellmap = HDGElasticity.CellMap([0.,0.],[1.,0.5])
coeffs = [1.0,-1.0,1.,1.]
update!(imap,coeffs)
normals = repeat([1.0,0.0],inner=(1,length(squad)))
HH = zeros(4,4)

facescale = HDGElasticity.scale_area(cellmap,normals)
HDGElasticity.HHop!(HH,sbasis,squad,imap,normals,facescale,2)
testHH = 0.25*testmatrix
testHH[2,:] .= 0.0
testHH[4,:] .= 0.0
@test allapprox(HH,testHH)

fill!(HH,0.0)
HDGElasticity.HHop!(HH,sbasis,squad,imap,normals,2facescale,2)
@test allapprox(HH,2testHH)

HH = HDGElasticity.HHop(sbasis,squad,imap,normals,normals,cellmap)
@test allapprox(HH,testHH)

ext: jl


text:
# Load Julia packages (libraries) needed for clip

using DrWatson
#@quickactivate "SR2StanPluto"
using StanSample
using StatisticalRethinking

df = CSV.read(sr_datadir("Howell1.csv"), DataFrame)

# Use only adults

df2 = filter(row -> row[:age] >= 18, df);
first(df2, 5)

heightsmodel = "
// Inferring the mean and std
data {
  int N;
  real<lower=0> h[N];
}
parameters {
  real<lower=0> sigma;
  real<lower=0,upper=250> mu;
}
model {
  // Priors for mu and sigma
  mu ~ normal(178, 0.1);
  sigma ~ uniform(0 , 50);

  // Observed heights
  h ~ normal(mu, sigma);
}
";

sm = SampleModel("heights", heightsmodel);

heightsdata = Dict("N" => length(df2[:, :height]), "h" => df2[:, :height]);

rc = stan_sample(sm, data=heightsdata);

if success(rc)
	println()
	q = read_samples(sm, :particles)
  q |> display
end

# End of `clip-31.jl`


text:
# ==============================================================================
#                                  BP model
# ==============================================================================
using SeisAcoustic

dir_work = joinpath(homedir(), "Desktop/BP");
dir_model= joinpath(dir_work, "physical_model");

dir_near = joinpath(dir_work, "near_boat");
dir_far  = joinpath(dir_work, "far_boat" );

path_vel  = joinpath(dir_model, "vel.rsf");
path_rho  = joinpath(dir_model, "rho.rsf");
path_src  = joinpath(dir_model, "wlet.rsf");

(hdr, vel)  = read_RSdata(path_vel);
(hdr, rho)  = read_RSdata(path_rho);
(hdr, wlet) = read_RSdata(path_src);

# reduced physical model
vel = vel[1:1500, 700:end-700];
rho = rho[1:1500, 700:end-700];
(nz, nx) = size(vel);
for i2 = 1 : nx
    for i1 = 1 : nz
        if vel[i1,i2] > 4510
           vel[i1,i2] = 4510
        end
    end
end

# vertical and horizontal grid size
dz = 6.25; dx = 6.25;

# time step size and maximum modelling length
dt = 0.0008; tmax = 14.0;

# towed streamer
nr = 481; delta_r = 12.5; near_offset = 150.0;
ns = 510; delta_s = 25.0;

# compute receiver and source index for SLO geometry
irz = Vector{Vector{Int64}}(undef, ns);
irx = Vector{Vector{Int64}}(undef, ns);
isz_near = Vector{Int64}(undef, ns);
isx_near = Vector{Int64}(undef, ns);
isz_far  = Vector{Int64}(undef, ns);
isx_far  = Vector{Int64}(undef, ns);
for i = 1 : ns
    l2 = (nx-1) * dx - (i-1)*delta_s
    l1 = l2 - (nr-1) * delta_r
    rec_loc = l1 : delta_r : l2
    irx[i] = floor.(Int64, rec_loc/dx) .+ 1
    irz[i] = 4 * ones(Int64, nr)
    isx_near[i] = floor(Int64, (l1-150     )/dx) + 1
    isz_near[i] = 3
    isx_far[i]  = floor(Int64, (l1-150-6000)/dx) + 1
    isz_far[i]  = 3
end

dir_obs = joinpath(dir_work, "near_boat");
get_shotgather_adaptive(dir_obs, isz_near, isx_near, 1e5*ricker(10.0,dt), irz, irx, vel, rho, dz, dx, dt, tmax,
                        location_flag="index", data_format=Float64, order=2, free_surface=true, npad=100)

# neurus last 100 (415-510)
i1 = 415; i2 = 510;
dir_obs = joinpath(dir_work, "far_boat");
get_shotgather_adaptive(dir_obs, isz_far[i1:i2], isx_far[i1:i2], 1e5*ricker(10.0,dt), irz[i1:i2], irx[i1:i2], vel, rho, dz, dx, dt, tmax,
                        location_flag="index", data_format=Float64, order=2, free_surface=true, npad=10)

# download data to local machine
scp -r wgao1@hijitus.physics.ualberta.ca:/home/wgao1/Desktop/BP/far_boat $HOME/Desktop/BP
scp -r wgao1@saig-ml.physics.ualberta.ca:/home/wgao1/Desktop/BP/near_boat $HOME/Desktop/synthetic_data/BP

# resample the recordings from 0.8 ms to 4ms
dir_work = joinpath(homedir(), "Desktop/BP");
dt_new = 0.004; dt = 0.0008;
interval = floor(Int64, dt_new/dt);

# near boat
ns = 510; # total number of shots
dir_near = joinpath(dir_work, "near_boat");
dir_near_resample = joinpath(dir_work, "near_boat_4ms");
mkdir(dir_near_resample)
for i =1 : ns
    println("processed $i shot")
    path_in  = join([dir_near "/recordings_" "$i" ".bin"]);
    path_out = join([dir_near_resample "/pressure_" "$i" ".rsf"]);
    rec      = read_recordings(path_in)
    p        = rec.p[1:interval:end,:]

    hdr      = RegularSampleHeader(p, d1=dt_new, d2=12.5)
    write_RSdata(path_out, hdr, p)
end

dir_near_resample = joinpath(dir_work, "near_boat_4ms");
path = joinpath(dir_near_resample, "pressure_1.rsf");
hdr  = read_RSheader(path);
p    = zeros(hdr.data_format, hdr.n1, hdr.n2, ns);
for i = 1 : ns
    path = join([dir_near_resample "/pressure_" "$i" ".rsf"])
    (hdr, d) = read_RSdata(path)
    p[:,:,i].= d
end
pathout = joinpath(dir_work, "near_boat_4ms.rsf");
hdr     = RegularSampleHeader(p, d1=dt_new, d2=12.5, d3=25.0, title="near boat shots");
write_RSdata(pathout, hdr, p);
rm(dir_near_resample, recursive=true);

# far boat
dir_far = joinpath(dir_work, "far_boat");
dir_far_resample = joinpath(dir_work, "far_boat_4ms");
mkdir(dir_far_resample)
for i = 1 : ns
    println("processed $i shot")
    path_in  = join([dir_far "/recordings_" "$i" ".bin"]);
    path_out = join([dir_far_resample "/pressure_" "$i" ".rsf"]);
    rec      = read_recordings(path_in)
    p        = rec.p[1:interval:end,:]

    hdr      = RegularSampleHeader(p, d1=dt_new, d2=12.5)
    write_RSdata(path_out, hdr, p)
end

# make a cube
path = joinpath(dir_far_resample, "pressure_1.rsf");
hdr  = read_RSheader(path);
p    = zeros(hdr.data_format, hdr.n1, hdr.n2, ns);
for i = 1 : ns
    path = join([dir_far_resample "/pressure_" "$i" ".rsf"])
    (hdr, d) = read_RSdata(path)
    p[:,:,i].= d
end
pathout = joinpath(dir_work, "far_boat_4ms.rsf");
hdr     = RegularSampleHeader(p, d1=dt_new, d2=12.5, d3=25.0, title="far boat shots");
write_RSdata(pathout, hdr, p);
rm(dir_far_resample, recursive=true);

# download data to local machine
scp -r wgao1@saig-ml.physics.ualberta.ca:/home/wgao1/Desktop/BP/near_boat_4ms.rsf $HOME/Desktop/synthetic_data/BP
scp -r wgao1@saig-ml.physics.ualberta.ca:/home/wgao1/Desktop/BP/far_boat_4ms.rsf $HOME/Desktop/synthetic_data/BP

# apply time shift to far boat data and formulate blended super-gather
dir_work   = joinpath(homedir(), "Desktop/synthetic_data/BP")
path_near  = joinpath(dir_work, "near_boat_4ms.rsf");
(hdr1, p1) = read_RSdata(path_near);

path_far   = joinpath(dir_work, "far_boat_4ms.rsf");
(hdr2, p2) = read_RSdata(path_far);

i3=1;
SeisPlotTX(p1[:,:,i3], cmap="gray", wbox=10, hbox=10);

i2 = 1;
SeisPlotTX(p1[:,i2,:], cmap="gray", wbox=10, hbox=10);

i3=1;
SeisPlotTX(p2[:,:,i3], cmap="gray", wbox=10, hbox=10);

i2 = 1;
SeisPlotTX(p2[:,i2,:], cmap="gray", wbox=10, hbox=10);

time_delay = -1.75 * rand(hdr1.n3)     # make distribution same as PGS data
hdr_time   = RegularSampleHeader(time_delay, title="time_delay applied to far boat");
path_time  = joinpath(dir_work, "time_delay.rsf");
write_RSdata(path_time, hdr_time, time_delay);

for i = 1 : hdr2.n3
    tmp = copy(p2[:,:,i])
    tmp1= time_shift(tmp, time_delay[i]*ones(hdr2.n2), hdr2.d1)
    p2[:,:,i] .= tmp1
    println("$i")
end
p1 .= p1 .+ p2;

i2 = 1;
SeisPlotTX(p1[:,i2,:], cmap="gray", wbox=10, hbox=10);

hdr_blend1 = RegularSampleHeader(p1, d1=hdr1.d1, d2=hdr1.d2, d3=hdr1.d3, title="blended data with near boat aligned");
path_blend1= joinpath(dir_work, "near_aligned.rsf");
write_RSdata(path_blend1, hdr_blend1, p1);

# make far boat data aligned
for i = 1 : hdr1.n2
    tmp = copy(p1[:,i,:])
    tmp1= time_shift(tmp, -time_delay, hdr1.d1)
    p1[:,i,:] .= tmp1
    println("$i")
end
i2 = 1;
SeisPlotTX(p1[:,i2,:], cmap="gray", wbox=10, hbox=10);

hdr_blend2 = RegularSampleHeader(p1, d1=hdr1.d1, d2=12.5, d3=25.0, title="blended data with far boat aligned");
path_blend2= joinpath(dir_work, "far_aligned.rsf");
write_RSdata(path_blend2, hdr_blend2, p1);

dir_pgs   = joinpath(homedir(), "Desktop/synthetic_data/PGS");
path_near = joinpath(dir_pgs, "data.rsf");
(hdr, d1) = read_RSdata(path_near);
dir_BP    = joinpath(homedir(), "Desktop/synthetic_data/BP");
path_near = joinpath(dir_BP, "near_aligned.rsf");
(hdr, p1) = read_RSdata(path_near);

i2 = 1;
SeisPlotTX(p1[:,i2,:], cmap="gray", wbox=10, hbox=10);
SeisPlotTX(d1[:,i2,:], cmap="gray", wbox=10, hbox=10);

# write them as a binary
dir_BP    = joinpath(homedir(), "Desktop/synthetic_data/BP");
path_near = joinpath(dir_BP, "near_aligned.rsf");
(hdr, p1_blend) = read_RSdata(path_near);

path_input_near = joinpath(dir_BP, "input_near_aligned.bin");
fid = open(path_input_near, "w");
write(fid, vec(p1_blend)); close(fid);

path_near = joinpath(dir_BP, "near_boat_4ms.rsf");
(hdr, p1_clear) = read_RSdata(path_near);

path_target_near = joinpath(dir_BP, "target_near_aligned.bin");
fid = open(path_target_near, "w");
write(fid, vec(p1_clear)); close(fid);

# write far boat data
dir_BP    = joinpath(homedir(), "Desktop/synthetic_data/BP");
path_far = joinpath(dir_BP, "far_aligned.rsf");
(hdr, p2_blend) = read_RSdata(path_far);

path_input_far = joinpath(dir_BP, "input_far_aligned.bin");
fid = open(path_input_far, "w");
write(fid, vec(p2_blend)); close(fid);

path_far = joinpath(dir_BP, "far_boat_4ms.rsf");
(hdr, p2_clear) = read_RSdata(path_far);

path_target_far = joinpath(dir_BP, "target_far_aligned.bin");
fid = open(path_target_far, "w");
write(fid, vec(p2_clear)); close(fid);

dir_BP     = joinpath(homedir(), "Desktop/synthetic_data/BP");
path_delay = joinpath(dir_BP, "time_delay.rsf");
(hdr, time_delay) = read_RSdata(path_delay);

path_time = joinpath(dir_BP, "time_delay.bin");
fid = open(path_time, "w");
write(fid, convert(vector{Float32}, vec(time_delay)); close(fid);

i2 = 1;
SeisPlotTX(p1_blend[:,i2,:], cmap="gray", wbox=10, hbox=10);
SeisPlotTX(p1_blend[:,i2,:]-p1_clear[:,i2,:], cmap="gray", wbox=10, hbox=10);

# ==============================================================================
#                                   Marmousi2
# ==============================================================================
# using SeisAcoustic
#
# dir_work = joinpath(homedir(), "Desktop/Marmousi2");
# dir_model= joinpath(dir_work, "physical_model_rsf");
# dir_near = joinpath(dir_work, "near_boat");
# dir_far  = joinpath(dir_work, "far_boat" );
#
# path_vp   = joinpath(dir_model, "vp.rsf");
# path_rho  = joinpath(dir_model, "rho.rsf");
# path_wlet = joinpath(dir_model, "wavelet_08ms.rsf");
#
# (hdr, vp)   = read_RSdata(path_vp);
# (hdr, rho)  = read_RSdata(path_rho);
# (hdr, wlet) = read_RSdata(path_wlet);
#
# # reduced physical model
# vp  = vp[1:2:2401, 1000:5:end];
# rho = rho[1:2:2401, 1000:5:end];
# (nz, nx) = size(vp);
#
# # top boundary condition
# free_surface = true;
#
# # vertical and horizontal grid size
# dz = 6.25; dx = 6.25;
#
# # time step size and maximum modelling length
# dt = hdr.d1; tmax = 16.0;
#
# # precision
# data_format = Float32;
# order       = 2;
#
# # organize these parameters into a structure
# fidiff = TdParams(rho, vp, free_surface, dz, dx, dt, tmax;
#                   data_format=data_format, order=order);
#
# # towed streamer
# nr = 257; delta_r = 12.5; near_offset = 150.0;
# ns = 257; delta_s = 25.0;
#
# # compute receiver and source index for SLO geometry
# irz = Vector{Vector{Int64}}(undef, ns);
# irx = Vector{Vector{Int64}}(undef, ns);
# isz_near = Vector{Int64}(undef, ns);
# isx_near = Vector{Int64}(undef, ns);
# isz_far  = Vector{Int64}(undef, ns);
# isx_far  = Vector{Int64}(undef, ns);
# for i = 1 : ns
#     l2 = (nx-1) * dx - (i-1)*delta_s
#     l1 = l2 - (nr-1) * delta_r
#     rec_loc = l1 : delta_r : l2
#     irx[i] = floor.(Int64, rec_loc/dx) .+ 1
#     irz[i] = 4 * ones(Int64, nr)
#     isx_near[i] = floor(Int64, (l1-150     )/dx) + 1
#     isz_near[i] = 3
#     isx_far[i]  = floor(Int64, (l1-150-6000)/dx) + 1
#     isz_far[i]  = 3
# end
#
# # # # near offset source
# srcs    = get_multi_sources(isz_near, isx_near, fidiff; p=wlet);
# dir_obs = joinpath(dir_work, "near_boat");
# get_observations(dir_obs, irz, irx, srcs, fidiff);
#
# # far offset source
# il = 206; iu=ns;
# srcs = get_multi_sources(isz_far[il:iu], isx_far[il:iu], fidiff; p=wlet);
# dir_obs = joinpath(dir_work, "far_boat");
# get_observations(dir_obs, irz[il:iu], irx[il:iu], srcs, fidiff);
#
# # download data to local machine
# scp -r wgao1@neurus.physics.ualberta.ca:/$HOME/Desktop/Marmousi2/far_boat $HOME/Desktop/Marmousi2/
# scp -r wgao1@saig-ml.physics.ualberta.ca:/home/wgao1/Desktop/Marmousi2/near_boat $HOME/Desktop/Marmousi2/
#
# # combine the file
# file_name = readdir("/Users/wenlei/Desktop/far_boat");
# dir_name  = joinpath(homedir(), "Desktop/rename")
# for i = 1 : length(file_name)
#     s = split(file_name[i], "_")[2]
#     a = parse(Int64, split(s, ".")[1]) + 205
#     path1 = joinpath("/Users/wenlei/Desktop/far_boat", file_name[i])
#     path2 = join([dir_name "/recordings_" "$a" ".bin"])
#     cp(path1, path2)
# end
#
# # resample the recordings from 0.8 ms to 4ms
# dir_work = joinpath(homedir(), "Desktop/Marmousi2");
# dt_new = 0.004; interval = floor(Int64, dt_new/dt)-1;
#
# dir_near = joinpath(dir_work, "near_boat");
# dir_near_4ms = joinpath(dir_work, "near_boat_4ms");
# for i = 1 : ns
#     path_in  = join([dir_near "/recordings_" "$i" ".bin"]);
#     path_out = join([dir_near_4ms "/pressure_" "$i" ".rsf"]);
#     rec      = read_recordings(path_in)
#     p        = rec.p[1:interval:end,:]
#
#     hdr      = RegularSampleHeader(p, d1=dt_new, d2=12.5)
#     write_RSdata(path_out, hdr, p)
# end
#
# dir_far = joinpath(dir_work, "far_boat");
# dir_far_4ms = joinpath(dir_work, "far_boat_4ms");
# for i = 1 : ns
#     path_in  = join([dir_far "/recordings_" "$i" ".bin"]);
#     path_out = join([dir_far_4ms "/pressure_" "$i" ".rsf"]);
#     rec      = read_recordings(path_in)
#     p        = rec.p[1:interval:end,:]
#
#     hdr      = RegularSampleHeader(p, d1=dt_new, d2=12.5)
#     write_RSdata(path_out, hdr, p)
# end
#
# # make a cube
# dir_work = joinpath(homedir(), "Desktop/Data_Beijing_Workshop/Marmousi2")
# dir_near_4ms = joinpath(dir_work, "near_boat_4ms");
# path = joinpath(dir_near_4ms, "pressure_1.rsf");
# hdr  = read_RSheader(path);
# p    = zeros(hdr.data_format, hdr.n1, hdr.n2, ns);
# for i = 1 : ns
#     path = join([dir_near_4ms "/pressure_" "$i" ".rsf"])
#     (hdr, d) = read_RSdata(path)
#     p[:,:,i].= d
# end
# path = joinpath(homedir(), "Desktop/Data_Beijing_Workshop/Marmousi2/near_4ms.rsf");
# hdr  = RegularSampleHeader(p, d1=0.004, d2=12.5, d3=25.0, unit1="s", unit2="meter", unit3="meter");
# write_RSdata(path, hdr, p);
#
# dir_far_4ms = joinpath(dir_work, "far_boat_4ms");
# path = joinpath(dir_far_4ms, "pressure_1.rsf");
# hdr  = read_RSheader(path);
# p    = zeros(hdr.data_format, hdr.n1, hdr.n2, ns);
# for i = 1 : ns
#     path = join([dir_far_4ms "/pressure_" "$i" ".rsf"])
#     (hdr, d) = read_RSdata(path)
#     p[:,:,i].= d
# end
# path = joinpath(homedir(), "Desktop/Marmousi2/far_4ms.rsf");
# hdr  = RegularSampleHeader(p, d1=0.004, d2=12.5, d3=25.0, unit1="s", unit2="meter", unit3="meter");
# write_RSdata(path, hdr, p);
#
#
# path1 = joinpath(homedir(), "Desktop/Marmousi2/near_4ms.rsf");
# path2 = joinpath(homedir(), "Desktop/Marmousi2/far_4ms.rsf");
# (hdr, p1) = read_RSdata(path1);
# (hdr, p2) = read_RSdata(path2);
#
# # manual blending the two cube
# delay_max = 1.75; tau = 1.75 * rand(ns);
# delay = zeros(Float64, ns);
# for i = 1 : ns
#     delay[i] = (round(Int64, tau[i]/dt)) * dt
# end
# hdr = RegularSampleHeader(delay);
# path = joinpath(homedir(), "Desktop/Marmousi2/delay_4ms.rsf");
# write_RSdata(path, hdr, delay);
#
# path = joinpath(homedir(), "Desktop/Marmousi2/delay_4ms.rsf");
# (hdr, delay) = read_RSdata(path);
# for i = 1 : ns
#     delay_idx = round(Int64, delay[i]/dt)
#     tmp = vcat(zeros(Float32, delay_idx, nr), p2[1:end-delay_idx,:,i])
#     p1[:,:,i] .= p1[:,:,i] .+ tmp
# end
# hdr = RegularSampleHeader(p1, d1=0.004, d2=12.5, d3=25.0, unit1="s", unit2="meter", unit3="meter");
# path = joinpath(homedir(), "Desktop/Marmousi2/blend_4ms.rsf");
# write_RSdata(path, hdr, p1);
#
#
# # save it as binary
# path_in = joinpath(homedir(), "Desktop/Marmousi2/near_4ms.rsf");
# (hdr, near) = read_RSdata(path_in);
# path_out = joinpath(homedir(), "Desktop/Marmousi2/near_4ms.bin");
# fid  = open(path_out, "w");
# write(fid, vec(near)); close(fid);
#
# path_in = joinpath(homedir(), "Desktop/Marmousi2/far_4ms.rsf");
# (hdr, far) = read_RSdata(path_in);
# path_out = joinpath(homedir(), "Desktop/Marmousi2/far_4ms.bin");
# fid  = open(path_out, "w");
# write(fid, vec(far)); close(fid);
#
# path_in = joinpath(homedir(), "Desktop/Marmousi2/blend_4ms.rsf");
# (hdr, blend) = read_RSdata(path_in);
# path_out = joinpath(homedir(), "Desktop/Marmousi2/blend_4ms.bin");
# fid  = open(path_out, "w");
# write(fid, vec(blend)); close(fid);
#
# path_in = joinpath(homedir(), "Desktop/Marmousi2/delay_4ms.rsf");
# (hdr, dither) = read_RSdata(path_in);
# path_out = joinpath(homedir(), "Desktop/Marmousi2/dither_4ms.bin");
# fid  = open(path_out, "w");
# write(fid, convert(Vector{Float32}, dither)); close(fid);
#
#
# # make 16ms sampling cube
# p1_down = p1[1:4:end,:,:];
# p2_down = p2[1:4:end,:,:];
# path1 = joinpath(homedir(), "Desktop/Marmousi2/near_16ms.rsf");
# path2 = joinpath(homedir(), "Desktop/Marmousi2/far_16ms.rsf");
# hdr1  = RegularSampleHeader(p1_down, d1=0.016, d2=12.5, d3=25.0, unit1="s", unit2="meter", unit3="meter");
# hdr2  = RegularSampleHeader(p2_down, d1=0.016, d2=12.5, d3=25.0, unit1="s", unit2="meter", unit3="meter");
# write_RSdata(path1, hdr1, p1_down);
# write_RSdata(path2, hdr2, p2_down);
#
# shot_idx = 1
# SeisPlotTX(p1[:,:,shot_idx], hbox=8, wbox=5, pclip=90, cmap="gray"); tight_layout();
# SeisPlotTX(p2[:,:,shot_idx], hbox=8, wbox=5, pclip=90, cmap="gray"); tight_layout();
#
# channel_idx = 20
# SeisPlotTX(p1[:,channel_idx,:], hbox=8, wbox=5, pclip=90, cmap="gray"); tight_layout();
# SeisPlotTX(p2[:,channel_idx,:], hbox=8, wbox=5, pclip=90, cmap="gray"); tight_layout();
#
# # manual blending the two cube
# delay_max = 1.75; tau = 1.75 * rand(ns);
# delay = zeros(Float64, ns);
# for i = 1 : ns
#     delay[i] = (round(Int64, tau[i]/dt)) * dt
# end
# hdr = RegularSampleHeader(delay);
# path = joinpath(homedir(), "Desktop/Marmousi2/delay_16ms.rsf");
# write_RSdata(path, hdr, delay);
#
# path = joinpath(homedir(), "Desktop/Marmousi2/delay_16ms.rsf");
# (hdr, delay) = read_RSdata(path);
# for i = 1 : ns
#     delay_idx = round(Int64, delay[i]/dt)
#     tmp = vcat(zeros(Float32, delay_idx, nr), p2[1:end-delay_idx,:,i])
#     p1[:,:,i] .= p1[:,:,i] .+ tmp
# end
# hdr = RegularSampleHeader(p1, d1=0.016, d2=12.5, d3=25.0, unit1="s", unit2="meter", unit3="meter");
# path = joinpath(homedir(), "Desktop/Marmousi2/blend_16ms.rsf");
# write_RSdata(path, hdr, p1);
#
# for i = 60 : 60 : ns
#     SeisPlotTX(p1[:,:,i], hbox=8, wbox=5, pclip=90, cmap="gray"); tight_layout();
# end
#
# for i = 40 : 40 : nr
#     SeisPlotTX(p1[:,i,:], hbox=8, wbox=5, pclip=90, cmap="gray"); tight_layout();
# end
#
#
# # examine the data by plotting
# idx = 10;
# path_rec = join([dir_near_4ms "/pressure_" "$idx" ".rsf"]);
# (hdr, p) = read_RSdata(path_rec);
# SeisPlotTX(p, hbox=8, wbox=5, pclip=90, cmap="gray",
#            dx=1, dy=0.004, ox=1, oy=0.0, xticks=50:50:nr, yticks=2:2:16, ticksize=15,
#            xlabel="Trace number", ylabel="Time (S)", labelsize=15); tight_layout();
#
# # amplitude spectra
# (f, a) = amplitude_spectra(p, 0.004);
# plot(f, a);
#
# iter = 200;
# a    = 0.95;
# x    = zeros(iter);
# for i = 1 : iter
#     x[i] = a^i
# end
# plot(x)
#
# rec = read_recordings(joinpath(dir_obs, "recordings_8.bin"))
# SeisPlotTX(rec.p)
#
# # # # near offset source
# srcs    = get_multi_sources(isz_near[1:2], isx_near[1:2], fidiff; p=wlet);
# dir_obs = joinpath(dir_work, "near_boat");
# get_shotgather(dir_obs, irz[1:2], irx[1:2], srcs, fidiff);
#
# # near offset source
# srcs = get_multi_sources(isz_far, isx_far, fidiff; p=wlet);
# dir_obs = joinpath(dir_work, "far_boat");
# get_observations(dir_obs, irz, irx, srcs, fidiff);


text:
#=##############################################################################
# # DESCRIPTION
#     Examples of VTK formatting and usage of VTKtools
# # AUTHORSHIP
#   * Author    : Eduardo J Alvarez
#   * Email     : Edo.AlvarezR@gmail.com
#   * Created   : Nov 2017
#   * License   : MIT License
=###############################################################################

function blade_example(; prompt=true, file_name="temp_blade00")
  airfoil_files = ["Cyl1.txt", "Cyl1.txt", "S815.txt", "S809.txt", "S826.txt"]
  # file_name = "temp_blade00"    # Name for vtk outpout file

  println("Defining variables...")
  prev_time = time()

  # PARAMETERS
  n_upper = 20             # Number of sections in upper surface of blade
  n_lower = 20          # Number of sections in lower surface of blade
  r = [1.0, 1.0, 8.0, 8.0, 8.0]           # Expansion ratio in both surfaces of each airfoil
  sections = [[(1.0, 10, 1.0, false)],      # Discretization between each airfoil
              [(1.0, 5, 1.0, false)],
              [(1.0, 15, 1.0, false)],
              [(1.0, 10, 1.0, false)],
      ]

  Rtip = 25.0           # Radius at blade tip
  Rhub = 1.0            # Radius of the hub
  pos = [0, 0.15, 0.2, 0.5, 1.0]           # Position along blade of each airfoil
  chords = [1.0, 0.6, 1.75, 3.0, 0.85]    # Chord length of each airfoil

  # Leading edge position of each airfoil
  Os = [[-chords[1]/2, 0, Rhub+pos[1]*(Rtip-Rhub)],
          [-chords[2]*7/8, 0, Rhub+pos[2]*(Rtip-Rhub)],
          [-chords[3]*4/6, 0, Rhub+pos[3]*(Rtip-Rhub)],
          [-chords[4]/4, 0, Rhub+pos[4]*(Rtip-Rhub)],
          [-chords[5]/8, 0, Rhub+pos[5]*(Rtip-Rhub)]
      ]
  # Orientation of chord of each airfoil (yaw, pitch, roll)
  orien = [ [0.0, 0.0, 0.0],
              [0.0, 0.0, 0.0],
              [10.0, 0.0, 0.0],
              [7.5, 0.0, 0.0],
              [0.0, 30.0, 0.0]
  ]

  crosssections = []        # It will store here the cross sections for lofting
  point_datas = []          # Dummy point data for good looking visuals

  println("\t Runtime: $(round(time()-prev_time,digits=1)) (s)")
  println("Processing airfoils...")
  prev_time = time()

  # Processes each airfoil geometry
  styles = ["--k", "--r", "--g", "--b", "--y", "--c"]
  org_points = []
  for (i,airfoil_file) in enumerate(airfoil_files)

      # Read airfoil file
      x,y = readcontour(data_path*airfoil_file; header_len=2)
      push!(org_points, [x,y])

      # Separate upper and lower sides to make the contour injective in x
      upper, lower = splitcontour(x, y)

      # Parameterize both sides independently
      fun_upper = gt.parameterize(upper[1], upper[2], zeros(Float64, size(upper[1])); inj_var=1)
      fun_lower = gt.parameterize(lower[1], lower[2], zeros(Float64, size(lower[1])); inj_var=1)

      # New discretization for both surfaces
      upper_points = gt.discretize(fun_upper, 0, 1, n_upper, r[1]; central=true)
      lower_points = gt.discretize(fun_lower, 0, 1, n_lower, r[1]; central=true)

      # Put both surfaces back together from TE over the top and from LE over the bottom.
      reverse!(upper_points)                           # Trailing edge over the top
      new_x = [point[1] for point in upper_points]
      new_y = [point[2] for point in upper_points]      # Leading edge over the bottom
      new_x = vcat(new_x, [point[1] for point in lower_points])
      new_y = vcat(new_y, [point[2] for point in lower_points])

      plot_airfoil(new_x, new_y; style=styles[i], label=airfoil_file)
  #     println("Close figure and press ENTER")
  #     readline()

      # Scales the airfoil acording to its chord length
      new_x = chords[i]*new_x
      new_y = chords[i]*new_y

      # Reformats into points
      npoints = size(new_x)[1]
      airfoil = Array{Float64, 1}[[new_x[j], new_y[j], 0] for j in 1:npoints]

      # Positions the airfoil along the blade in the right orientation
      Oaxis = gt.rotation_matrix(orien[i][1], orien[i][2], orien[i][3])
      invOaxis = inv(Oaxis)
      airfoil = gt.countertransform(airfoil, invOaxis, Os[i])

      push!(crosssections, airfoil)
      push!(point_datas, [j for j in npoints*(i-1).+1:npoints*i])
  end

  println("\t Runtime: $(round(time()-prev_time,digits=1)) (s)")
  println("Lofting and generating VTK files...")
  prev_time = time()

  # Generates cells in VTK Legacy format
  out = gt.multilines2vtkmulticells(crosssections, sections;
                                        point_datas=point_datas)
  points, vtk_cells, point_data = out

  # Formats the point data for generateVTK
  data = []
  push!(data, Dict(
              "field_name" => "Point_index",
              "field_type" => "scalar",
              "field_data" => point_data
              )
   )

  # Generates the vtk file
  gt.generateVTK(file_name, points; cells=vtk_cells, point_data=data)

  println("\t Runtime: $(round(time()-prev_time,digits=1)) (s)")

  # Calls paraview
  run(`paraview --data="$(file_name).vtk;"`)

  # Deletes files
  if prompt
   print("Delete vtk files? ([y]/n) ")
   y = readline()
  else
   y = "y"
  end
  if y=="y"; run(`rm -f $(file_name).vtk`); end;
end


text:
#=

A chat application using both HttpServer and HTTP to do
the same thing: Start a new task for each browser (tab) that connects.

To use:
    - include("chat_explore.jl") in REPL
    - start a browser on address 127.0.0.1:8000, and another on 127.0.0.1:8080
    - inspect global variables starting with 'last' while the chat is running asyncronously

To call in from other devices, figure out your IP address on the network and change the 'gatekeeper' code.

Note that type of 'lastreq' changes depending on whether the last call was made through HttpServer or HTTP.

Functions used as arguments are explicitly defined with names instead of anonymous functions (do..end constructs).
This may improve debugging readability at the cost of increased verbosity.

=#
global lastreq = 0
global lastws= 0
global lastmsg= 0
global lastws= 0

using HttpServer
using HTTP
using WebSockets
const CLOSEAFTER = Base.Dates.Second(1800)
const HTTPPORT = 8080
const HTTPSERVERPORT = 8000
const URL = "127.0.0.1"
const USERNAMES = Dict{String, WebSocket}()
const HTMLSTRING = readstring(Pkg.dir("WebSockets","examples","chat_explore.html"));

# If we are to access a websocket from outside
# it's websocket handler function, we need some kind of
# mutable container for storing references:
const WEBSOCKETS = Dict{WebSocket, Int}()

"""
Called by 'gatekeeper', this function stays active while the
particular websocket is open. The argument is an open websocket.
Other instances of the function run in other tasks. The tasks
are generated by either HTTP or HttpServer.
"""
function coroutine(thisws)
    global lastws = thisws
    push!(WEBSOCKETS, thisws => length(WEBSOCKETS) +1 )
    t1 = now() + CLOSEAFTER
    username = ""
    while now() < t1
        # This next call waits for a message to
        # appear on the socket. If there is none,
        # this task yields to other tasks.
        data, success = readguarded(thisws)
        !success && break
        global lastmsg = msg = String(data)
        print("Received: $msg ")
        if username == ""
            username = approvedusername(msg, thisws)
            if username != ""
                println("from new user $username ")
                !writeguarded(thisws, username) && break
                println("Tell everybody about $username")
                foreach(keys(WEBSOCKETS)) do ws
                    writeguarded(ws, username * " enters chat")
                end
            else
                println(", username taken!")
                !writeguarded(thisws, "Username taken!") && break
            end
        else
            println("from $username ")
            distributemsg(msg, thisws)
            startswith(msg, "exit") && break
        end
    end
    exitmsg = username == "" ? "unknown" : username * " has left"
    distributemsg(exitmsg, thisws)
    println(exitmsg)
    # No need to close the websocket. Just clean up external references:
    removereferences(thisws)
    nothing
end

function removereferences(ws)
    haskey(WEBSOCKETS, ws) && pop!(WEBSOCKETS, ws)
    for (discardname, wsref) in USERNAMES
        if wsref === ws
            pop!(USERNAMES, discardname)
            break
        end
    end
    nothing
end

function approvedusername(msg, ws)
    !startswith(msg, "userName:") && return ""
    newname = msg[length("userName:") + 1:end]
    newname =="" && return ""
    haskey(USERNAMES, newname) && return ""
    push!(USERNAMES, newname => ws)
    newname
end

function distributemsg(msgout, not_to_ws)
    foreach(keys(WEBSOCKETS)) do ws
        if ws !== not_to_ws
            writeguarded(ws, msgout)
        end
    end
    nothing
end

"""
`Server => gatekeeper(Request, WebSocket) => coroutine(WebSocket)`

The gatekeeper makes it a little harder to connect with
malicious code. It inspects the request that was upgraded
to a a websocket.
"""
function gatekeeper(req, ws)
    global lastreq = req
    global lastws = ws
    orig = WebSockets.origin(req)
    if startswith(orig, "http://localhost") || startswith(orig, "http://127.0.0.1")
        coroutine(ws)
    else
        warn("Unauthorized websocket connection, $orig not approved by gatekeeper")
    end
    nothing
end

"Request to response. Response is the predefined HTML page with some javascript"
req2resp(req::HttpServer.Request, resp) = HTMLSTRING |> Response
req2resp(req::HTTP.Request) =       HTMLSTRING |> HTTP.Response

# Both server definitions need two function wrappers; one handler function for page requests,
# one for opening websockets (which the javascript in the HTML page will try to do)
server_httpserver = Server(HttpHandler(req2resp), WebSocketHandler(gatekeeper))
server_HTTP = WebSockets.ServerWS(HTTP.HandlerFunction(req2resp), WebSockets.WebsocketHandler(gatekeeper))

# Start the HTTP server asyncronously, and stop it later
litas_HTTP = @schedule WebSockets.serve(server_HTTP, URL, HTTPPORT)
@schedule begin
    println("HTTP server listening on $URL:$HTTPPORT for $CLOSEAFTER")
    sleep(CLOSEAFTER.value)
    println("Time out, closing down $HTTPPORT")
    Base.throwto(litas_HTTP, InterruptException())
end

# Start the HttpServer asyncronously, stop it later
litas_httpserver = @schedule run(server_httpserver, HTTPSERVERPORT)
@schedule begin
          println("HttpServer listening on $URL:$HTTPSERVERPORT for $CLOSEAFTER")
          sleep(CLOSEAFTER.value + 2)
          println("Time out, closing down $HTTPSERVERPORT")
          Base.throwto(litas_httpserver, InterruptException())
end

nothing


text:
using CairoMakie, SpecialFunctions
p(s, β) = (β + 1)*a(β) * s^β * exp(-a(β)*s^(β+1))
a(β) = gamma(((β + 2)/(β + 1))^(β + 1))
s = LinRange(0,3,30)
colors = ["#FF410D", "#6EE2FF", "#F7C530", "#95CC5E", "#D0DFE6", "#F79D1E"]

with_theme(theme_dark()) do
    fig = Figure(resolution = (600,400))
    ax = Axis(fig[1,1]; palette = (; color = colors))
    for β in [0,1], step in [:pre, :center, :post]
        stairs!(s, p.(s, β), step = step, label = "$(β), :$(step)")
    end
    lines!(s, p.(s, 1), color = :grey90)
    lines!(s, p.(s, 0), color = :grey70)
    text!(L"p(s)=(\beta+1)a_{\beta}\,s^{\beta}\exp(-a_{\beta}s^{\beta+1})",
        position = (0.7, 0.95), color = :white)
    text!(L"a_{\beta} = \Gamma[(\beta+2)/(\beta+1)]^{\beta + 1}",
        position = (1.7, 0.25), color = :white)
    axislegend("β, step", position = :rt)
    save("stairs.png", fig, px_per_unit = 2)
    display(fig)
end


text:
function TriangularDistribution(arg0::jdouble, arg1::jdouble, arg2::jdouble)
    return TriangularDistribution((jdouble, jdouble, jdouble), arg0, arg1, arg2)
end

function cumulative_probability(obj::TriangularDistribution, arg0::jdouble)
    return jcall(obj, "cumulativeProbability", jdouble, (jdouble,), arg0)
end

function density(obj::TriangularDistribution, arg0::jdouble)
    return jcall(obj, "density", jdouble, (jdouble,), arg0)
end

function get_mode(obj::TriangularDistribution)
    return jcall(obj, "getMode", jdouble, ())
end

function get_numerical_mean(obj::TriangularDistribution)
    return jcall(obj, "getNumericalMean", jdouble, ())
end

function get_numerical_variance(obj::TriangularDistribution)
    return jcall(obj, "getNumericalVariance", jdouble, ())
end

function get_support_lower_bound(obj::TriangularDistribution)
    return jcall(obj, "getSupportLowerBound", jdouble, ())
end

function get_support_upper_bound(obj::TriangularDistribution)
    return jcall(obj, "getSupportUpperBound", jdouble, ())
end

function inverse_cumulative_probability(obj::TriangularDistribution, arg0::jdouble)
    return jcall(obj, "inverseCumulativeProbability", jdouble, (jdouble,), arg0)
end

function is_support_connected(obj::TriangularDistribution)
    return jcall(obj, "isSupportConnected", jboolean, ())
end


text:
#=
Copyright (c) 2018-2022 Chris Coey, Lea Kapelevich, and contributors

This Julia package Hypatia.jl is released under the MIT license; see LICENSE
file in the root directory or at https://github.com/chriscoey/Hypatia.jl

sets of native test instances
=#

inst_preproc = [
    "dimension1",
    "consistent1",
    "inconsistent1",
    "inconsistent2",
    ]

inst_infeas = [
    "primalinfeas1",
    "primalinfeas2",
    "primalinfeas3",
    "dualinfeas1",
    "dualinfeas2",
    "dualinfeas3",
    ]

inst_minimal = [
    "primalinfeas3",
    "dualinfeas3",
    "epinorminf4",
    "hyporootdettri4",
    ]

inst_cones_few = [
    "nonnegative1",
    "possemideftri1",
    "possemideftri5",
    "doublynonnegativetri1",
    "possemideftrisparse2",
    "possemideftrisparse5",
    "linmatrixineq1",
    "epinorminf4",
    "epinorminf6",
    "epinormeucl1",
    "epipersquare1",
    "epinormspectral1",
    "matrixepipersquare2",
    "generalizedpower1",
    "hypopowermean1",
    "hypogeomean1",
    "hyporootdettri1",
    "hypoperlog1",
    "hypoperlogdettri1",
    "epipersepspectral_vector1",
    "epipersepspectral_matrix2",
    "epirelentropy4",
    "epitrrelentropytri1",
    "wsosinterpnonnegative1",
    "wsosinterpnonnegative4",
    "wsosinterppossemideftri1",
    "wsosinterpepinormone1",
    "wsosinterpepinormeucl2",
    ]

# superset of inst_cones_few
inst_cones_many = [
    "nonnegative1",
    "nonnegative2",
    "nonnegative3",
    "nonnegative4",
    "possemideftri1",
    "possemideftri2",
    "possemideftri3",
    "possemideftri4",
    "possemideftri5",
    "possemideftri6",
    "possemideftri7",
    "possemideftri8",
    "possemideftri9",
    "doublynonnegativetri1",
    "doublynonnegativetri2",
    "doublynonnegativetri3",
    "possemideftrisparse1",
    "possemideftrisparse2",
    "possemideftrisparse3",
    "possemideftrisparse4",
    "possemideftrisparse5",
    "linmatrixineq1",
    "linmatrixineq2",
    "linmatrixineq3",
    "epinorminf1",
    "epinorminf2",
    "epinorminf3",
    "epinorminf4",
    "epinorminf5",
    "epinorminf6",
    "epinorminf7",
    "epinormeucl1",
    "epinormeucl2",
    "epinormeucl3",
    "epipersquare1",
    "epipersquare2",
    "epipersquare3",
    "epipersquare4",
    "epinormspectral1",
    "epinormspectral2",
    "epinormspectral3",
    "epinormspectral4",
    "matrixepipersquare1",
    "matrixepipersquare2",
    "matrixepipersquare3",
    "generalizedpower1",
    "generalizedpower2",
    "generalizedpower3",
    "generalizedpower4",
    "hypopowermean1",
    "hypopowermean2",
    "hypopowermean3",
    "hypopowermean4",
    "hypopowermean5",
    "hypopowermean6",
    "hypogeomean1",
    "hypogeomean2",
    "hypogeomean3",
    "hypogeomean4",
    "hypogeomean5",
    "hypogeomean6",
    "hyporootdettri1",
    "hyporootdettri2",
    "hyporootdettri3",
    "hyporootdettri4",
    "hypoperlog1",
    "hypoperlog2",
    "hypoperlog3",
    "hypoperlog4",
    "hypoperlog5",
    "hypoperlog6",
    "hypoperlog7",
    "hypoperlogdettri1",
    "hypoperlogdettri2",
    "hypoperlogdettri3",
    "hypoperlogdettri4",
    "epipersepspectral_vector1",
    "epipersepspectral_vector2",
    "epipersepspectral_vector3",
    "epipersepspectral_vector4",
    "epipersepspectral_matrix1",
    "epipersepspectral_matrix2",
    "epipersepspectral_matrix3",
    "epirelentropy1",
    "epirelentropy2",
    "epirelentropy3",
    "epirelentropy4",
    "epirelentropy5",
    "epitrrelentropytri1",
    "epitrrelentropytri2",
    "epitrrelentropytri3",
    "epitrrelentropytri4",
    "wsosinterpnonnegative1",
    "wsosinterpnonnegative2",
    "wsosinterpnonnegative3",
    "wsosinterpnonnegative4",
    "wsosinterpnonnegative5",
    "wsosinterppossemideftri1",
    "wsosinterppossemideftri2",
    "wsosinterppossemideftri3",
    "wsosinterpepinormone1",
    "wsosinterpepinormone2",
    "wsosinterpepinormone3",
    "wsosinterpepinormeucl1",
    "wsosinterpepinormeucl2",
    "wsosinterpepinormeucl3",
    ]

inst_indirect = [
    "indirect1",
    "indirect2",
    "indirect3",
    "indirect4",
    "indirect5",
    ]


text:
# Autogenerated wrapper script for qt5_jll for x86_64-apple-darwin14
export libqt5quick, libqt5quickcontrols2, libqt5widgets, libqt5svg, libqt5core

using Xorg_libX11_jll
using Xorg_libXext_jll
using Xorg_glproto_jll
using Xorg_libxcb_jll
using Xorg_xcb_util_wm_jll
using Xorg_xcb_util_image_jll
using Xorg_xcb_util_keysyms_jll
using Xorg_xcb_util_renderutil_jll
using xkbcommon_jll
using Libglvnd_jll
using Fontconfig_jll
using Glib_jll
## Global variables
PATH = ""
LIBPATH = ""
LIBPATH_env = "DYLD_FALLBACK_LIBRARY_PATH"

# Relative path to `libqt5quick`
const libqt5quick_splitpath = ["lib", "QtQuick.framework", "Versions", "5", "QtQuick"]

# This will be filled out by __init__() for all products, as it must be done at runtime
libqt5quick_path = ""

# libqt5quick-specific global declaration
# This will be filled out by __init__()
libqt5quick_handle = C_NULL

# This must be `const` so that we can use it with `ccall()`
const libqt5quick = "@rpath/QtQuick"

# Relative path to `libqt5quickcontrols2`
const libqt5quickcontrols2_splitpath = ["lib", "QtQuickControls2.framework", "Versions", "5", "QtQuickControls2"]

# This will be filled out by __init__() for all products, as it must be done at runtime
libqt5quickcontrols2_path = ""

# libqt5quickcontrols2-specific global declaration
# This will be filled out by __init__()
libqt5quickcontrols2_handle = C_NULL

# This must be `const` so that we can use it with `ccall()`
const libqt5quickcontrols2 = "@rpath/QtQuickControls2"

# Relative path to `libqt5widgets`
const libqt5widgets_splitpath = ["lib", "QtWidgets.framework", "Versions", "5", "QtWidgets"]

# This will be filled out by __init__() for all products, as it must be done at runtime
libqt5widgets_path = ""

# libqt5widgets-specific global declaration
# This will be filled out by __init__()
libqt5widgets_handle = C_NULL

# This must be `const` so that we can use it with `ccall()`
const libqt5widgets = "@rpath/QtWidgets"

# Relative path to `libqt5svg`
const libqt5svg_splitpath = ["lib", "QtSvg.framework", "Versions", "5", "QtSvg"]

# This will be filled out by __init__() for all products, as it must be done at runtime
libqt5svg_path = ""

# libqt5svg-specific global declaration
# This will be filled out by __init__()
libqt5svg_handle = C_NULL

# This must be `const` so that we can use it with `ccall()`
const libqt5svg = "@rpath/QtSvg"

# Relative path to `libqt5core`
const libqt5core_splitpath = ["lib", "QtCore.framework", "Versions", "5", "QtCore"]

# This will be filled out by __init__() for all products, as it must be done at runtime
libqt5core_path = ""

# libqt5core-specific global declaration
# This will be filled out by __init__()
libqt5core_handle = C_NULL

# This must be `const` so that we can use it with `ccall()`
const libqt5core = "@rpath/QtCore"

"""
Open all libraries
"""
function __init__()
    global artifact_dir = abspath(artifact"qt5")

    # Initialize PATH and LIBPATH environment variable listings
    global PATH_list, LIBPATH_list
    # We first need to add to LIBPATH_list the libraries provided by Julia
    append!(LIBPATH_list, [joinpath(Sys.BINDIR, Base.LIBDIR, "julia"), joinpath(Sys.BINDIR, Base.LIBDIR)])
    # From the list of our dependencies, generate a tuple of all the PATH and LIBPATH lists,
    # then append them to our own.
    foreach(p -> append!(PATH_list, p), (Xorg_libX11_jll.PATH_list, Xorg_libXext_jll.PATH_list, Xorg_glproto_jll.PATH_list, Xorg_libxcb_jll.PATH_list, Xorg_xcb_util_wm_jll.PATH_list, Xorg_xcb_util_image_jll.PATH_list, Xorg_xcb_util_keysyms_jll.PATH_list, Xorg_xcb_util_renderutil_jll.PATH_list, xkbcommon_jll.PATH_list, Libglvnd_jll.PATH_list, Fontconfig_jll.PATH_list, Glib_jll.PATH_list,))
    foreach(p -> append!(LIBPATH_list, p), (Xorg_libX11_jll.LIBPATH_list, Xorg_libXext_jll.LIBPATH_list, Xorg_glproto_jll.LIBPATH_list, Xorg_libxcb_jll.LIBPATH_list, Xorg_xcb_util_wm_jll.LIBPATH_list, Xorg_xcb_util_image_jll.LIBPATH_list, Xorg_xcb_util_keysyms_jll.LIBPATH_list, Xorg_xcb_util_renderutil_jll.LIBPATH_list, xkbcommon_jll.LIBPATH_list, Libglvnd_jll.LIBPATH_list, Fontconfig_jll.LIBPATH_list, Glib_jll.LIBPATH_list,))

    global libqt5quick_path = normpath(joinpath(artifact_dir, libqt5quick_splitpath...))

    # Manually `dlopen()` this right now so that future invocations
    # of `ccall` with its `SONAME` will find this path immediately.
    global libqt5quick_handle = dlopen(libqt5quick_path)
    push!(LIBPATH_list, dirname(libqt5quick_path))

    global libqt5quickcontrols2_path = normpath(joinpath(artifact_dir, libqt5quickcontrols2_splitpath...))

    # Manually `dlopen()` this right now so that future invocations
    # of `ccall` with its `SONAME` will find this path immediately.
    global libqt5quickcontrols2_handle = dlopen(libqt5quickcontrols2_path)
    push!(LIBPATH_list, dirname(libqt5quickcontrols2_path))

    global libqt5widgets_path = normpath(joinpath(artifact_dir, libqt5widgets_splitpath...))

    # Manually `dlopen()` this right now so that future invocations
    # of `ccall` with its `SONAME` will find this path immediately.
    global libqt5widgets_handle = dlopen(libqt5widgets_path)
    push!(LIBPATH_list, dirname(libqt5widgets_path))

    global libqt5svg_path = normpath(joinpath(artifact_dir, libqt5svg_splitpath...))

    # Manually `dlopen()` this right now so that future invocations
    # of `ccall` with its `SONAME` will find this path immediately.
    global libqt5svg_handle = dlopen(libqt5svg_path)
    push!(LIBPATH_list, dirname(libqt5svg_path))

    global libqt5core_path = normpath(joinpath(artifact_dir, libqt5core_splitpath...))

    # Manually `dlopen()` this right now so that future invocations
    # of `ccall` with its `SONAME` will find this path immediately.
    global libqt5core_handle = dlopen(libqt5core_path)
    push!(LIBPATH_list, dirname(libqt5core_path))

    # Filter out duplicate and empty entries in our PATH and LIBPATH entries
    filter!(!isempty, unique!(PATH_list))
    filter!(!isempty, unique!(LIBPATH_list))
    global PATH = join(PATH_list, ':')
    global LIBPATH = join(LIBPATH_list, ':')

    # Add each element of LIBPATH to our DL_LOAD_PATH (necessary on platforms
    # that don't honor our "already opened" trick)
    #for lp in LIBPATH_list
    #    push!(DL_LOAD_PATH, lp)
    #end
end  # __init__()


text:
#=
Description
    Module for pairwise UAV encounters. The dynamics implemented are
    based on the simple kinematic model, as seen in the UAVs module.
    You can modify the discretization scheme by varying XDIM, YDIM,
    PDIM, and VDIM. The current values are chosen for a reasonable
    compute time for quick visualization purposes.
=#

module DoubleUAVs

using POMDPs

export NSTATES, COC, DoubleUAV, get_next_state, gen_reward

const DT = 5.0              # [s]
const DTI = 1.0             # [s]
const G = 9.8               # [m/s^2]

const MIN_SEP_SQ = 500.0^2  # [m]
const RBF_STD_DIST = 500.0  # [m]
const RBF_INV_VAR = 1 / RBF_STD_DIST^2   # [1/m^2]

const PEN_ACTION = 0.02
const PEN_CLOSENESS = 10
const PEN_MIN_SEP = 1000
const PEN_CONFLICT = 1.0

const STATE_DIM = 5
const ACTION_DIM = 2

const X = 1                 # [m] relative x-position
const Y = 2                 # [m] relative y-position
const P = 3                 # [rad] relative heading
const V1 = 4                # [m/s] ac1 speed
const V2 = 5                # [m/s] ac2 speed

const XDIM = 11
const XMIN = -2e3
const XMAX = 2e3

const YDIM = 11
const YMIN = -2e3
const YMAX = 2e3

const PDIM = 5
const PMIN = 0.0
const PMAX = 2 * pi

const VDIM = 3
const VMIN = 10
const VMAX = 20

xs = linspace(XMIN, XMAX, XDIM)
ys = linspace(YMIN, YMAX, YDIM)
ps = linspace(PMIN, PMAX, PDIM)
vs = linspace(VMIN, VMAX, VDIM)

const NSTATES = XDIM * YDIM * PDIM * VDIM^2 + 1
tmp_states = zeros(STATE_DIM, NSTATES)
istate = 1
for iv2 = 1:VDIM
    for iv1 = 1:VDIM
        for ip = 1:PDIM
            for iy = 1:YDIM
                for ix = 1:XDIM
                    tmp_states[X, istate] = xs[ix]
                    tmp_states[Y, istate] = ys[iy]
                    tmp_states[P, istate] = ps[ip]
                    tmp_states[V1, istate] = vs[iv1]
                    tmp_states[V2, istate] = vs[iv2]
                    istate += 1
                end # for ix
            end # for iy
        end # for ip
    end # for iv1
end # for iv2

const TERM_STATE_VAR = 1e5
const TERM_STATE = TERM_STATE_VAR * ones(STATE_DIM)
tmp_states[:, end] = TERM_STATE
const STATES_OBSERVATIONS = tmp_states

const COC = -1
indivActions = deg2rad([-20, -10, 0, 10, 20, rad2deg(COC)])
nIndivActions = length(indivActions)
actions = zeros(2, nIndivActions^2)
iaction = 1
for ia2 = 1:nIndivActions
    for ia1 = 1:nIndivActions
        actions[1, iaction] = indivActions[ia1]
        actions[2, iaction] = indivActions[ia2]
        iaction += 1
    end # for ia1
end # for ia2
const ACTIONS = actions

const SIGMA_V = 2.0                 # [m/s]
const SIGMA_B = deg2rad(4.0)        # [rad]
const SIGMA_B_COC = deg2rad(10.0)   # [rad]
const SIGMAS = [0 SIGMA_B -SIGMA_B 0 0 0 0 0 0;
                0 0 0 SIGMA_B -SIGMA_B 0 0 0 0;
                0 0 0 0 0 SIGMA_V -SIGMA_V 0 0;
                0 0 0 0 0 0 0 SIGMA_V -SIGMA_V]
const WEIGHTS = [1 / 3, 2 / (3 * (size(SIGMAS, 2) - 1)) *
                             ones(size(SIGMAS, 2) - 1)]
const SIGMA_DEF = zeros(size(SIGMAS, 1))

type DoubleUAV

    pomdp::     POMDP
    sigmas::    Matrix{Float64}
    weights::   Vector{Float64}

    function DoubleUAV()
        p = POMDP(STATES_OBSERVATIONS, ACTIONS, STATES_OBSERVATIONS, TERM_STATE)
        return new(p, SIGMAS, WEIGHTS)
    end # function DoubleUAV

end # DoubleUAV

function norm_angle(angle::Float64)
    return ((angle % (2 * pi)) + 2 * pi) % (2 * pi)
end # function norm_angle

function gen_next_state(dt::Float64 = DT)
    function get_next_state(state::Vector{Float64}, iaction::Int64, isigma::Int64=1)
        if state[1] == TERM_STATE_VAR
            return TERM_STATE
        else
            x = state[1]
            y = state[2]
            p2 = state[3]
            v1 = state[4]
            v2 = state[5]

            b1 = ACTIONS[1, iaction]
            b2 = ACTIONS[2, iaction]

            # for sigma-point sampling
            if b1 != COC
                b1 = b1 + SIGMAS[1, isigma]
            else
                b1 = sign(SIGMAS[1, isigma]) * SIGMA_B_COC
            end # if

            if b2 != COC
                b2 = b2 + SIGMAS[2, isigma]
            else
                b2 = sign(SIGMAS[2, isigma]) * SIGMA_B_COC
            end # if

            v1 = v1 + SIGMAS[3, isigma]
            v2 = v2 + SIGMAS[4, isigma]

            if b1 == 0 || b2 == 0 # straight-line path(s)
                if b2 != 0 # uav1 straight path
                    dp2 = dt * G * tan(b2) / v2
                    r2 = abs(v2^2 / (G * tan(b2)))
                    xr = x + r2 * sign(dp2) * (sin(p2) - sin(p2 - dp2)) - v1 * dt
                    yr = y + r2 * sign(dp2) * (-cos(p2) + cos(p2 - dp2))
                    pr = norm_angle(p2 + dp2)

                    if xr < XMIN || xr > XMAX ||
                       yr < YMIN || yr > YMAX
                        return TERM_STATE
                    else
                        return [xr, yr, pr, v1, v2]
                    end # if
                elseif b1 != 0 # uav2 straight path
                    dp1 = dt * G * tan(b1) / v1
                    r1 = abs(v1^2 / (G * tan(b1)))
                    x = x + v2 * dt * cos(p2) - r1 * sign(dp1) * sin(dp1)
                    y = y + v2 * dt * sin(p2) - r1 * sign(dp1) * (cos(dp1) - 1)
                    xr = x * cos(dp1) + y * sin(dp1)
                    yr = -x * sin(dp1) + y * cos(dp1)
                    pr = norm_angle(p2 - dp1)

                    if xr < XMIN || xr > XMAX ||
                       yr < YMIN || yr > YMAX
                        return TERM_STATE
                    else
                        return [xr, yr, pr, v1, v2]
                    end # if
                else # both straight paths
                    xr = x + v2 * dt * cos(p2) - v1 * dt
                    yr = y + v2 * dt * sin(p2)
                    pr = norm_angle(p2)

                    if xr < XMIN || xr > XMAX ||
                       yr < YMIN || yr > YMAX
                        return TERM_STATE
                    else
                        return [xr, yr, pr, v1, v2]
                    end # if
                end # if
            else # both curved paths
                dp1 = dt * G * tan(b1) / v1
                dp2 = dt * G * tan(b2) / v2
                r1 = abs(v1^2 / (G * tan(b1)))
                r2 = abs(v2^2 / (G * tan(b2)))
                x = x + r2 * sign(dp2) * (sin(p2) - sin(p2 - dp2)) - r1 * sign(dp1) * sin(dp1)
                y = y + r2 * sign(dp2) * (-cos(p2) + cos(p2 - dp2)) - r1 * sign(dp1) * (-1 + cos(dp1))
                xr = x * cos(dp1) + y * sin(dp1)
                yr = -x * sin(dp1) + y * cos(dp1)
                pr = norm_angle(p2 + dp2 - dp1)

                if xr < XMIN || xr > XMAX ||
                   yr < YMIN || yr > YMAX
                    return TERM_STATE
                else
                    return [xr, yr, pr, v1, v2]
                end # if
            end # if
        end # if
    end # function get_next_state
    return get_next_state
end # function gen_next_state

function gen_reward(pen_conflict::Float64=PEN_CONFLICT)
    get_next_state = gen_next_state(DTI)
    function get_reward(state::Vector{Float64}, iaction::Int64)
        action = ACTIONS[:, iaction]

        reward = 0

        if action[1] == COC
            action[1] = 0.0
        else
            reward = -pen_conflict
        end # if

        if action[2] == COC
            action[2] = 0.0
        else
            reward = reward - pen_conflict
        end # if

        if state[1] == TERM_STATE_VAR
            return reward - PEN_ACTION * rad2deg(norm(action))^2
        else
            reward = reward - PEN_ACTION * rad2deg(norm(action))^2

            minSepSq = Inf
            for ti = 1:DT / DTI
                minSepSq = min(minSepSq, norm(state)^2)
                state = get_next_state(state, iaction)
            end # for ti

            if minSepSq < MIN_SEP_SQ
                reward = reward - PEN_MIN_SEP
            end # if

            reward = reward - PEN_CLOSENESS * exp(-minSepSq * RBF_INV_VAR)
            return reward
        end # if
    end # function get_reward
    return get_reward
end # function gen_reward

get_next_state = gen_next_state(DT)

end # module DoubleUAVs


text:
struct CachedDict
    cachefile::String
    d
    CachedDict(fn::Function, filename::String) = begin
        new(filename, fromcacheorcreate(fn, filename))
    end
end
Base.haskey(c::CachedDict, key) = haskey(c.d, key)
Base.getindex(c::CachedDict, key) = c.d[key]
Base.setindex!(c::CachedDict, value, key) = update!(c, value, key, true)

function cachedict(c::CachedDict)
    if !isdir(dirname(c.cachefile))
        mkpath(dirname(c.cachefile))
    end
    serialize(c.cachefile, c.d)
end

function update!(c::CachedDict, value, key, writetodisc = false)
    c.d[key] = value
    writetodisc && cachedict(c)
    value
end

function fromcacheorcreate(fn, filename)
    !isfile(filename) && return fn()
    try
        deserialize(filename)
    catch _err
        return fn()
    end
end


text:
function correlate(w)

 # the model correlations for quantities F, B, S, dedf based on del and E

  del = w[:,1]
  E = w[:,2]./del .- 1

  F = 4.8274*E.^4 - 5.9816*E.^3 + 4.0274*E.^2 + 0.23247.*E .+ 0.15174

  B = 0.5*(-225.86.*E.^3 - 3016.6.*E.^2 - 208.68*E .- 17.915
                + 131.9*E.^3 - 167.32*E.^2 + 76.642.*E .- 11.068)
  B[E.<-0.0616] = -225.86.*E[E.<-0.0616].^3 - 3016.6.*E[E.<-0.0616].^2 - 208.68*E[E.<-0.0616] .- 17.915
  B[E.>-0.0395] =  131.9*E[E.>-0.0395].^3 - 167.32*E[E.>-0.0395].^2 + 76.642.*E[E.>-0.0395] .- 11.068

  S = 0.5*(451.55*E.^3 + 2010*E.^2 + 138.96*E .+ 11.296
                        - 96.739*E.^3 + 117.74*E.^2 - 46.432*E .+ 6.8074)
  S[E.<-0.0582] = 451.55*E[E.<-0.0582].^3 + 2010*E[E.< -0.0582].^2 + 138.96*E[E.< -0.0582] .+ 11.296
  S[E.>-0.042]  = -96.739*E[E.>-0.042].^3 + 117.74*E[E.> -0.042].^2 - 46.432*E[E.> -0.042] .+ 6.8074

  dfde = 4*4.8274*E.^3 - 3*5.9816*E.^2 + 2*4.0274*E .+ 0.23247

  return del, E, F ,B, S, dfde

end

max_stars_repo_name: desanga/IBMFVM

max_issues_repo_name: desanga/IBMFVM

max_forks_repo_name: desanga/IBMFVM


text:
using Blink
using Base.Test

cleanup = !AtomShell.isinstalled()

cleanup && AtomShell.install()

# open window and wait for it to initialize
w = Window(Blink.@d(:show => false)); sleep(5.0)

# make sure the window is really active
@test_approx_eq @js(w, Math.log(10)) log(10)

@test sprint(Blink.jsexpr, :(Dict("a" => 1, :b => 10))) == "{\"a\":1,b:10}"

cleanup && AtomShell.remove()


text:
using BenchmarkTools, Expokit, ExpmV

BenchmarkTools.DEFAULT_PARAMETERS.seconds = 120

const p = 1e-3

bench = Dict{String,Any}[]
for d in 2.^(5:10)
    println("d = $d ...")
    b1 = @benchmarkable expm(full_r)*rv setup=((full_r,rv) = (full(randn()*sprandn($d,$d,p/2)+1im*sprandn($d,$d,p/2)), rv=normalize(randn($d)+1im*randn($d))))
    b2 = @benchmarkable Expokit.expmv(rt,r,rv) setup=((rt,r,rv)=(rand(), sprandn($d,$d,p/2)+1im*sprandn($d,$d,p/2), normalize(randn($d)+1im*randn($d))))
    b3 = @benchmarkable ExpmV.expmv(rt,r,rv) setup=((rt,r,rv)=(rand(), sprandn($d,$d,p/2)+1im*sprandn($d,$d,p/2), normalize(randn($d)+1im*randn($d))))

    t1 = run(b1)
    t2 = run(b2)
    t3 = run(b3)

    push!(bench,Dict("d"=>d, "p"=>p,"expm"=>t1,"expokit"=>t2,"expmv"=>t3))
end

JLD.save("master-bench.jld","bench",bench)


text:
function test_random_normal_layer(backend::Backend, T, eps)
  println("-- Testing RandomNormal Layer on $(typeof(backend)){$T}...")

  ############################################################
  # Prepare Data for Testing
  ############################################################

  batch_sizes = [100,200]
  N = length(batch_sizes)
  tensor_dim = abs(rand(Int)) % 4 + 1
  output_dims = abs(rand(Int, tensor_dim)) % 3 + 1
  println("    > Random output ", tensor_dim, output_dims)

  ############################################################
  # Setup
  ############################################################

  # batch size is determined by
  layer = RandomNormalLayer(tops = [:data1, :data2], output_dims=output_dims,
                            eltype=T, batch_sizes=batch_sizes)
  state = setup(backend, layer, Blob[], Blob[])

  layer_data = [Array(T, tuple(output_dims..., batch_sizes[i]))
                for i in 1:N]

  forward(backend, state, Blob[])
  for i in 1:N
      copy!(layer_data[i], state.blobs[i])
      @test (abs(mean(layer_data[i])) < 4e-1)
      @test all(-1000 .< layer_data[i] .< 1000)
  end

    # we should have sample from zero mean, unit stddev in state.blobs[1]
    # the size should be as expected

    # output should be different on subsequent calls
  layer_data2 = [Array(T, tuple(output_dims..., batch_sizes[i]))
                for i in 1:N]

  forward(backend, state, Blob[])
  for i in 1:N
      copy!(layer_data2[i], state.blobs[i])
      @test (abs(mean(layer_data2[i])) < 4e-1)
      @test all(-1000 .< layer_data2[i] .< 1000)
      @test 0 < norm(vec(layer_data[i] - layer_data2[i])) < 5sqrt(prod(output_dims)*batch_sizes[i])
  end

  shutdown(backend, state)
end
function test_random_normal_layer(backend::Backend)
  test_random_normal_layer(backend, Float32, 1e-5)
  test_random_normal_layer(backend, Float64, 1e-10)
end

if test_gpu
  test_random_normal_layer(backend_gpu)
end

if test_cpu
  test_random_normal_layer(backend_cpu)
end


text:
using HDF5,JLD
using Docile
using Lumberjack

import JLD.JldFile
# if for whatever reason the backend needs to change
# then adding write/delete/update methods here with a
# different backend / file_handle type would be added

abstract StoreBackend

# possibly allows for filehandle in here and read/write privilidge?
# this is to avoid closing and opening handles all the time for
# many individual writes? bad idea? probably so don't put in for now.
# writes should generally be done more batch like - everything modified
# in memory as much as possible - then written to a track
# We want to try different forms of hdf5 e.g. compressed and mmapped, these can't be mixed

type JldStoreBackend <: StoreBackend
    path::AbstractString
    description::AbstractString
    compressed::Bool
    mmapped::Bool
end

# hdf5-blosc, hdf5-gzip hdf5-mmap, hdf5-uncompressed (uncompressed)
function getdb(dbpath::AbstractString; dbtype = "hdf5", description="" )
    # check for suffix
    if dbtype == "hdf5" # default with compression
        return JldStoreBackend(dbpath,description,true,false)
    elseif dbtype == "hdf5-blosc"
        return JldStoreBackend(dbpath,description,true,false)
    elseif dbtype == "hdf5-gzip"
        return JldStoreBackend(dbpath,description,true,false)
    elseif dbtype == "hdf5-mmap"
        return JldStoreBackend(dbpath,description,false,true)
    elseif dbtype == "hdf5-uncompressed"
       return JldStoreBackend(dbpath,description,false,false)
    end
end

"""
Sets up filehandle and writes values to a genomic track

# Arguments
   * `file_handle::JldBackend`: backend store type
   * `seq_id::String`: identifying name of genomic sequence
   * `track_id::String`: identifying name of track
   * `values::Number`
   * `overwrite::Boolean` : keyword argument indicating whether to overwrite the track

# Return
  * JldBackend
"""
function write_track{T <: Number}(backend_store::JldStoreBackend,
                                  track_id::AbstractString, seq_id::AbstractString,
                                  values::Vector{T}; overwrite=false )
    genomic_store_path = backend_store.path
    fid=nothing
    if !isfile(genomic_store_path)
        fid=jldopen(genomic_store_path,"w",compress=backend_store.compressed)
    else
        fid=jldopen(genomic_store_path,"r+",compress=backend_store.compressed)
    end

    if !has(fid.plain,"$seq_id/$track_id") & !overwrite
        Lumberjack.info("Writing $seq_id")
        _write_track(fid, track_id, seq_id,values)
        Lumberjack.info("Finished writing $seq_id")
    else
        Lumberjack.info("Track: $seq_id/$track_id exists ... skipping")
    end
    close(fid)
end

function delete_track!(backend_store::JldStoreBackend,track_id::AbstractString, seq_id::AbstractString)
    genomic_store_path = backend_store.path
    fid=jldopen(genomic_store_path,"r+",compress=backend_store.compressed)
    Lumberjack.info("delete $seq_id")
    _delete_track!(fid, track_id, seq_id)
    Lumberjack.info("Finished removing $seq_id")
    close(fid)
end

function read_track(backend_store::JldStoreBackend,
                                 track_id::AbstractString,  seq_id::AbstractString,
                                 start_pos::Integer, stop_pos::Integer)
    genomic_store_path = backend_store.path
    fid=nothing
    try
        fid=h5open(genomic_store_path,"r")
    catch e
        Lumberjack.error("$e\n Trying to open $genomic_store_path")
    end
    if has(fid,"$seq_id/$track_id")
        if start_pos == 0 && stop_pos == 0
            track_slice=h5read(genomic_store_path,bytestring("$seq_id/$track_id"))
        else
            track_slice=h5read(genomic_store_path,bytestring("$seq_id/$track_id"),(start_pos:stop_pos,))
        end
        close(fid)
        return track_slice
    end
    close(fid)
    return []
end

function update_track!{T <: Number}(backend_store::JldStoreBackend,
                                  track_id::AbstractString,  seq_id::AbstractString,
                                  start_pos::AbstractString, stop_pos::AbstractString,
                                  values::Vector{T})
    Lumberjack.info("Updating $seq_id")
    genomic_store_path = backend_store.path
    fid=jldopen(genomic_store_path,"r+",compress=backend_store.compressed)
    _udpate_track!(fid,track_id , seq_id , start_pos,stop_pos,values)
    Lumberjack.info("Finished updating $seq_id")
    close(fid)
end

function _istrack(backend_store::JldStoreBackend,track_id::AbstractString)
    genomic_store_path = backend_store.path
    fid=nothing
    try
        fid=h5open(genomic_store_path,"r")
    catch e
        Lumberjack.error("$e\n Trying to open $genomic_store_path")
    end
    for n in names(fid)
        if has(fid,"$n/$track_id")
            close(fid)
            return true
        end
    end
    close(fid)
    return false
end

"""
  _hastrack( backend_store::JldStoreBackend,track_id::AbstractString )

  A track is a set of intervals or points associated with a sequence.
  So we look through all sequence names for evidence that a track exists.
  Should replace _istrack() _istrack will be deprecated in favour of this
  method.
"""

function _hastrack( backend_store::JldStoreBackend,track_id::AbstractString  )
    # _hastrack
    _istrack(backend_store,track_id)
end

"""
   _hasseqtrack(  backend_store::JldStoreBackend,track_id::AbstractString, seq_id::AbstractString )
   A track is a set of intervals or points associated with a sequence.
   This checks the existance of a seq_id/track_id path.
"""

function _hasseqtrack(  backend_store::JldStoreBackend,track_id::AbstractString, seq_id::AbstractString )
    genomic_store_path = backend_store.path
    fid=nothing
    try
         fid=h5open(genomic_store_path,"r")
    catch e
        Lumberjack.error("$e\n Trying to open $genomic_store_path")
    end
    if has(fid,"$seq_id/$track_id")
       close(fid)
       return true
    end
    close(fid)
    return false
end

# backend store - hdf5
# call these _write_track_hdf5  (crud-hdf5.jl?)
# catch errors here to allow logging of errors
# with specific messages relevant to these functions

function _write_track{T <: Number}(file_handle::JldFile, track_id::AbstractString, seq_id::AbstractString, values::Vector{T} )
    try
        write(file_handle, bytestring("$seq_id/$track_id"), values)
    catch err
        Lumberjack.error("Problem writing track: $err")
    end
end

function _delete_track!(file_handle::JldFile,track_id::AbstractString, seq_id::AbstractString)
    try
        delete!(file_handle,bytestring("$seq_id/$track_id"))
    catch err
        Lumberjack.error("Problem removing track: $err")
    end
end

function _update_track!{T <: Number}(file_handle::JldFile, track_id::AbstractString, seq_id::AbstractString,
                                      start_pos::Integer,stop_pos::Integer,values::Vector{T} )
    dset = file_handle["$seq_id/$track_id"]
    if stop_pos - start_pos != length(values)
        error("coordinates do not match length of values")
    end
    try
        dset[start_pos:stop_pos] = values
    catch err
        Lumberjack.error("Problem updating track: $err")
    end
end

"""
    _list_tracks

    Data organised as:

    "chr1"
     -> TRACKID1
     -> TRACKID2

    "chr2"
      -> TRACKID1
      -> TRACKID2

     Check whether a track_id exists under a sequence id.

"""
function _list_sequence_ids(backend_store::JldStoreBackend,track_id::AbstractString)
    genomic_store_path = backend_store.path
    fid=nothing
    try
        fid=h5open(genomic_store_path,"r")
    catch e
        Lumberjack.error("$e\n Trying to open $genomic_store_path")
    end
    seq_ids=[]
    for n in names(fid)
        if has(fid,"$n/$track_id")
            push!(seq_ids,n)
        end
    end
   close(fid)
   return seq_ids
end

function _list_tracks(backend_store::JldStoreBackend)
end

# a search function - e.g. indices or values  > than some value or matching a value .. prob not at this level - better done somewhere else
# specialise could be a missing value thing
# missing value tracks are very similar to storing intervals as 1 or 0
# this should all be done at a higher level than crud.


text:
"""
	vertCongruence(V::Lar.Array{Float64,2}; ϵ=1e-6)::Tuple{Lar.Points, Array{Array{Int,1},1}}

Evaluates the Vertex Congruence for 3D-points ``V ∈ ℳ(3,n)``.

The function determines the points of `V` closer than `ϵ` and
builds a new Vertex Set made of the representative of each point cluster.

The method returns:
 - the new Vertex Set
 - a map that, for every new vertex, identifies the old vertices it is made of
"""
function vertCongruence(V; ϵ=1e-6)
	Vcls    = []
	let visited = [],  kdtree = NearestNeighbors.KDTree(V)
		for vidx = 1 : size(V, 2)  if !(vidx in visited)
			nearvs = NearestNeighbors.inrange(kdtree, V[:, vidx], ϵ)
			push!(Vcls, nearvs)
			push!(visited, nearvs...)
		end  end
	end

	V = hcat([sum(V[:, cl], dims=2)/length(cl) for cl in Vcls]...)
	return V, Vcls
end


text:
using ProblemParser
using Utils

using Chain
using Lazy
using Underscores

@Base.kwdef struct Point
    x::Int = 0
    y::Int = 0
end

struct Line
    a::Point
    b::Point
end

horizontal(line::Line) = line.a.x == line.b.x
vertical(line::Line) = line.a.y == line.b.y
axisaligned(line) = horizontal(line) || vertical(line)

const file = find_input(@__FILE__)
content = @> file slurp
width = height = 0
const input = per_line(content, false) do line
    x1, y1, x2, y2 = @_ line |>
        split(__, " -> ") |>
        map(split(_, ","), __) |>
        Iterators.flatten |>
        map(curry(parse, Int), __)
    global width = max(width, x1, x2)
    global height = max(height, y1, y2)
    Line(Point(x1,y1), Point(x2,y2))
end
height += 1
width += 1

function sirange(start, stop)
    start > stop && return range(start, stop, step=-1)
    range(start, stop)
end

function drawAA(input)
    field = zeros(Int, height, width)
    for line in input
        same = line.a
        if horizontal(line)
            field[sirange(line.a.y, line.b.y), same.x] .+= 1
        elseif vertical(line)
            field[same.y, sirange(line.a.x, line.b.x)] .+= 1
        end
    end
    field
end

function drawdiag!(field, input)
    for line in filter(!axisaligned, input)
        for (x,y) in zip(sirange(line.a.x, line.b.x), sirange(line.a.y, line.b.y))
            field[y, x] += 1
        end
    end
end

field = drawAA(input)
answer1 = count(>=(2), field)
assertequal(answer1, 8350)

drawdiag!(field, input)
answer2 = count(>=(2), field)
assertequal(answer2, 19374)

max_stars_repo_name: CmdQ/AoC


text:
abstract type Galaxy end

struct Spiral <: Galaxy
	name::Symbol
	diskHeight::Number
	radius::Number
	units::Quantity
end


text:
using SerializingChannels
using Test
using TypeSerializers

@testset "SerializingChannels" begin
    integer_value = 42
    c = Channel{Vector{UInt8}}(1)
    integer_channel = SerializingChannel{Int, Channel{Vector{UInt8}}, SerializationTypeSerializer{Int}}(c)
    @test !isready(integer_channel)
    @test !isready(c)
    put!(integer_channel, integer_value)
    @test isready(integer_channel)
    @test isready(c)
    @test take!(integer_channel) == integer_value
end


text:
function EstimateAsymptoticVarianceSampleMean(X::AbstractArray{Float64,1}, alpha::Float64 = 1/3)
  # estimate the asymptotic variance for a stationary time series using block means
  # see v.d. Vaart, Time Series, Section 5.1
  # alpha determines block size through l = [n^alpha]
  n = length(X)
  l = Int(floor(n^alpha))
  Z = [mean(X[i:(i+l-1)]) for i in 1:(n-l+1)]
  mu = mean(X)
  v = l/(n-l+1) * sumabs2(Z-mu)
  F(x) = 1/(n-l+1) * sum(sqrt(l)*(Z-mu).<=x)
  return (v,F)
end

function BatchMeans(X::AbstractArray{Float64,1}; alpha::Float64 = 1/3, n_batches = -1)

  n = length(X)
  if (n_batches < 0)
    batch_size = Int(floor(n^alpha))
    n_batches = Int(floor(n / n^alpha))
  else
    batch_size = Int(floor(n/n_batches))
  end
  if n_batches > n/2
    return Inf
  end
  unused = n - batch_size * n_batches
  batch_means = [mean(X[unused + (i-1)*batch_size + 1:unused + i*batch_size]) for i in 1:n_batches]
  v = var(batch_means)*batch_size
end

function ess(X::AbstractArray{Float64,1};alpha::Float64=1/3, variance = -1, n_batches=-1)

  asvar = BatchMeans(X,alpha=alpha,n_batches=n_batches)
  if variance < 0
    variance = var(X)
  end
  return length(X) * variance/asvar
end

function SquaredJumpingDistance(X::AbstractArray{Float64})

  if ndims(X) == 1
    n_samples = length(X)
    differences = X[2:end] - X[1:end-1]
  else
    n_samples = size(X)[2]
    differences = X[:,2:end] - X[:,1:end-1]
  end
  return sumabs2(differences)/(n_samples-1)

end

function PooledVariance(AsVar::Vector{Float64}, n_samples::Vector{Int}, SampleVariance::Vector{Float64})

  return sum(n_samples .* SampleVariance ./ AsVar)/sum(n_samples ./AsVar)

end

max_stars_repo_name: jbierkens/zigzag-experiments

max_issues_repo_name: jbierkens/zigzag-experiments

max_forks_repo_name: jbierkens/zigzag-experiments


text:
using Test
using FlightMechanics

ψ, θ, ϕ = 0.5, 0.05, 0.2  # rad
p, q, r = 1.0, 0.2, 0.1  # rad/s

# Test regression
ear = body_angular_velocity_to_euler_angles_rates(p, q, r, θ, ϕ)
pqr = euler_angles_rates_to_body_angular_velocity(ear..., θ, ϕ)
@test isapprox(pqr, [p, q, r])

# Test against trimmer function
ψ_dot = 0.2
pqr1 = euler_angles_rates_to_body_angular_velocity(ψ_dot, 0.0, 0.0, θ, ϕ)
pqr2 = turn_rate_angular_velocity(ψ_dot, θ, ϕ)
@test isapprox(pqr1, pqr2)


text:
using ReinforcementLearningBase
using ReinforcementLearningCore
using Random
using Test
using StatsBase
using Distributions: probs, Normal, logpdf
using ReinforcementLearningEnvironments
using Flux
using Zygote
using CUDA

@testset "ReinforcementLearningCore.jl" begin
    include("core/core.jl")
    include("components/components.jl")
    include("utils/utils.jl")
    include("extensions.jl")
end


text:
# ## Animation
z = get_maximal_state(storage)
z = [[z[1] for t = 1:100]..., z..., [z[end] for t = 1:100]...]
T = length(z)
anim = MeshCat.Animation(convert(Int, floor(1.0 / 0.01)))
build_robot(mech, vis=vis)
for t = 1:T
    MeshCat.atframe(anim, t) do
        set_robot(vis, mech, z[t])
    end
end
MeshCat.setanimation!(vis, anim)
set_camera!(vis, cam_pos=[3,-0,0], zoom=0.84)

set_floor!(vis, x=0.2, y=10, z=0.02, color=RGBA(0,0,0,1))

# ## HMTL scene
render_static(vis)
open(joinpath(@__DIR__, "atlas_drop.html"), "w") do file
    write(file, static_html(vis))
end


text:
"""
    Sound

Module that exports wrappers around audio methods in `PortAudio`.
* `sound` and `soundsc` : methods for audio playback.
* `record` : method for audio recording.
"""
module Sound

include("soundsc.jl")
include("record.jl")

# support SampledSignals iff user has loaded that package
function __init__()
    @require SampledSignals = "bd7594eb-a658-542f-9e75-4c4d8908c167" include("sample-buf.jl")
end

# Phase vocoder
include("hann.jl")
include("findpeaks.jl")
include("phase_vocoder.jl")

end # module


text:
export fftshift_view, ifftshift_view
export fftshift2d, ifftshift2d, fftshift2d_view, ifftshift2d_view

# add constructor here, might be merged to ShiftedArrays
# this prevents that CircShiftedArrays get nested with twice application
# https://github.com/JuliaArrays/ShiftedArrays.jl/pull/44
function ShiftedArrays.CircShiftedArray(csa::CircShiftedArray, n = Tuple(0 for i in 1:N))
    CircShiftedArray(parent(csa), n .+ csa.shifts)
end

"""
    fftshift_view(A [, dims])

Result is semantically equivalent to `fftshift(A, dims)` but returns
a view instead.
"""
function fftshift_view(mat::AbstractArray{T, N}, dims=ntuple(identity, Val(N))) where {T, N}
    ShiftedArrays.circshift(mat, ft_center_diff(size(mat), dims))
end

"""
    ifftshift_view(A [, dims])

Result is semantically equivalent to `fftshift(A, dims)` but returns
a view instead.
"""
function ifftshift_view(mat::AbstractArray{T, N}, dims=ntuple(identity, Val(N))) where {T, N}
    diff = .-(ft_center_diff(size(mat), dims))
    return ShiftedArrays.circshift(mat, diff)
end

"""
    rfftshift_view(A, dims)

Shifts the frequencies to the center expect for `dims[1]` because there os no negative
and positive frequency.
"""
function rfftshift_view(mat::AbstractArray{T, N}, dims=ntuple(identity, Val(N))) where {T, N}
    ShiftedArrays.circshift(mat, rft_center_diff(size(mat), dims))
end

"""
    irfftshift_view(A, dims)
Shifts the frequencies back to the corner except for `dims[1]` because there os no negative
and positive frequency.
"""
function irfftshift_view(mat::AbstractArray{T, N}, dims=ntuple(identity, Val(N))) where {T, N}
    ShiftedArrays.circshift(mat ,.-(rft_center_diff(size(mat), dims)))
end

"""
    fftshift2d(mat::AbstractArray{T, N}) where {T, N}
    Short-hand for fftshift(mat, (1,2)) performing only a 2D ft.
    See fft for details.
"""
function fftshift2d(mat::AbstractArray{T, N}) where {T, N}
    fftshift(mat,(1,2))
end

"""
    ifftshift2d(mat::AbstractArray{T, N}) where {T, N}
    Short-hand for ifftshift(mat, (1,2)) performing only a 2D inverse ft.
    See ifft for details.
"""
function ifftshift2d(mat::AbstractArray{T, N}) where {T, N}
    ifftshift(mat,(1,2))
end

"""
    fftshift2d_view(mat::AbstractArray{T, N}) where {T, N}
    Short-hand for fftshift_view(mat, (1,2)) performing only a 2D ft.
    See fft for details.
"""
function fftshift2d_view(mat::AbstractArray{T, N}) where {T, N}
    fftshift_view(mat,(1,2))
end

"""
    ifftshift2d_view(mat::AbstractArray{T, N}) where {T, N}
    Short-hand for ifftshift_view(mat, (1,2)) performing only a 2D inverse ft.
    See ifft for details.
"""
function ifftshift2d_view(mat::AbstractArray{T, N}) where {T, N}
    ifftshift_view(mat,(1,2))
end


text:
"""
    ellipsoid_pointwise(T, r, shift)

Compute all lattice points in the ellipsoid ``\\{ n∈ \\mathbb{Z}^g \\,:\\, ||v(n)|| < r, v(n) = √π T(n+[[Y⁻¹y]])\\}``.

# Arguments
- `T::Array{<:Real}`: basis of lattice
- `r::Real`: radius of ellipsoid
- `shift::Array{<:Real}`: ``[[Y⁻¹y]]``.
"""
function ellipsoid_pointwise(T::Array{<:Real}, r::Real, shift::Array{<:Real})
    L = Lattice(T);
    z = zeros(Int64, L.m);
    points = enum_ellipsoid(L, z, L.m, r^2/π, 0.0, shift);
    if ~(z in points)
        push!(points, z);
    end
    return points;
end

"""
    ellipsoid_uniform(T, r)

Compute all lattice points in the deformed ellipsoid ``\\{ n ∈\\mathbb{Z}^g\\,:\\, π (n-c)^t Y (n-c) < r^2\\,,
|c\\_j|<1\\,,\\forall j=1,\\ldots,g\\}``.
"""
function ellipsoid_uniform(T::Array{<:Real}, r::Real)
    L = Lattice(T);
    # compute union of pointwise ellipsoids centered at the vertices of a unit hypercube of dimension m
    pointset = Set();
    for i = 0:2^L.m-1
        shift = [(-1)^j for j in digits(i, base = 2, pad = L.m)];
        points = enum_ellipsoid(L, zeros(Int64,L.m), L.m, r^2/π, 0.0, shift);
        union!(pointset, points);
    end
    union!(pointset, [zeros(Int64, L.m)]);
    return collect(pointset);
end

"""
    enum_ellipsoid(L, x, i, radius, l, shift)

Enumerates all lattice points in subtree at height i, where parent nodes are defined by x. Helper method for [`ellipsoid_pointwise`](@ref) and [`ellipsoid_uniform`](@ref).
"""
function enum_ellipsoid(L::Lattice, x::Array{<:Integer}, i::Integer, radius::Real, l::Real, shift::Array{<:Real})
    points = Array{Int64}[];
    sol = copy(x);
    c = -sum(Float64[(x[j]+shift[j])*L.μ[j,i] for j=i+1:L.m]) - shift[i];
    interval = sqrt((radius-l)/L.r[i]); # interval of possibilities for x[i]
    for k = floor(Int64, c - interval) : ceil(Int64, c + interval)
        l_new = l + (k-c)^2*L.r[i];
        if l_new <= radius # possible new solution found
            sol[i] = k;
            if i == 1 & ~(l_new == 0) # leaf of enumeration tree
                push!(points, copy(sol));
            elseif i > 1 # solve subtree enumeration recursively
                subtree_points = enum_ellipsoid(L, sol, i-1, radius, l_new, shift);
                append!(points, subtree_points);
            end
        end
    end
    return points;
end


text:
immutable VDPInitDist end
sampletype(::Type{VDPInitDist}) = TagState
function rand(rng::AbstractRNG, d::VDPInitDist)
    return TagState([0.0, 0.0], 8.0*rand(rng, 2)-4.0)
end

initial_state_distribution(::VDPTagProblem) = VDPInitDist()


text:
import JSON2

include("client.jl")

```
    Resource(kind::String,
    version::String,
    group::String,
    namespaced::Bool,
    schema::Dict{String, Any})

A Resource represents a Kubernetes resource, it can be a CRD or core resource.
https://kubernetes.io/docs/concepts/extend-kubernetes/api-extension/custom-resources/
```
mutable struct Resource
    kind::String
    plural::String
    version::String
    group::String
    namespaced::Bool
    schema::Dict{String, Any}
    type
end

function parse(filepath::String)::Resource
end

function install(client::K8sClient, resource::Resource)
end

function get(client::K8sClient, resource::Resource, name::String)::Dict{String, Any}
    uri = ""
    if resource.group == "core"
        uri = joinpath("/api/", resource.version, resource.plural, name)
    else
        uri = joinpath("/apis/", resource.group, resource.version, resource.plural, name)
    end
    resp = request(client, "GET", uri)
    str = String(resp.body)
    jobj = JSON2.Parser.parse(str)
    return jobj
end

function get(client::K8sClient, resource::Resource, name::String, namespace::String)::Dict{String, Any}
    uri = ""
    if resource.group == "core"
        uri = joinpath("/api/", resource.version, "namespaces", namespace, resource.plural, name)
    else
        uri = joinpath("/apis/", resource.group, resource.version, "namespaces", namespace, resource.plural, name)
    end
    resp = request(client, "GET", uri)
    str = String(resp.body)
    jobj = JSON2.Parser.parse(str)
    return jobj
end

function list(client::K8sClient, resource::Resource; labels::Dict{String, Any}=Nothing)::Dict{String, Any}
    uri = ""
    if resource.group == "core"
        uri = joinpath("/api/", resource.version, resource.plural)
    else
        uri = joinpath("/apis/", resource.group, resource.version, resource.plural)
    end
    resp = request(client, "GET", uri)
    str = String(resp.body)
    jobj = JSON2.Parser.parse(str)
    return jobj
end

function list(client::K8sClient, resource::Resource, namespace::String; labels::Dict{String, Any}=Nothing)::Dict{String, Any}
    uri = ""
    if resource.group == "core"
        uri = joinpath("/api/", resource.version, "namespaces", namespace, resource.plural)
    else
        uri = joinpath("/apis/", resource.group, resource.version, "namespaces", namespace, resource.plural)
    end
    resp = request(client, "GET", uri)
    str = String(resp.body)
    jobj = JSON2.Parser.parse(str)
    return jobj
end

function create(client::K8sClient, resource::Resource, namespace::String, spec::Dict{String, Any})::Dict{String, Any}
    uri = ""
    if resource.group == "core"
        uri = joinpath("/api/", resource.version, "namespaces", namespace, resource.plural)
    else
        uri = joinpath("/apis/", resource.group, resource.version, "namespaces", namespace, resource.plural)
    end
    bod = JSON2.write(spec)
    resp = request(client, "POST", uri, bod)
    str = String(resp.body)
    jobj = JSON2.Parser.parse(str)
    return jobj
end

function create(client::K8sClient, resource::Resource, spec::Dict{String, Any})::Dict{String, Any}
    uri = ""
    if resource.group == "core"
        uri = joinpath("/api/", resource.version, resource.plural)
    else
        uri = joinpath("/apis/", resource.group, resource.version, resource.plural)
    end
    bod = JSON2.write(spec)
    resp = request(client, "POST", uri, bod)
    str = String(resp.body)
    jobj = JSON2.Parser.parse(str)
    return jobj
end

function update(client::K8sClient, resource::Resource, spec::Dict{String, Any})::Dict{String, Any}
    uri = ""
    if resource.group == "core"
        uri = joinpath("/api/", resource.version, resource.plural)
    else
        uri = joinpath("/apis/", resource.group, resource.version, resource.plural)
    end
    bod = JSON2.write(spec)
    resp = request(client, "PATCH", uri, bod)
    str = String(resp.body)
    jobj = JSON2.Parser.parse(str)
    return jobj
end

function update(client::K8sClient, resource::Resource, namespace::String, spec::Dict{String, Any})::Dict{String, Any}
    uri = ""
    if resource.group == "core"
        uri = joinpath("/api/", resource.version, "namespaces", namespace, resource.plural)
    else
        uri = joinpath("/apis/", resource.group, resource.version, "namespaces", namespace, resource.plural)
    end
    bod = JSON2.write(spec)
    resp = request(client, "PATCH", uri, bod)
    str = String(resp.body)
    jobj = JSON2.Parser.parse(str)
    return jobj
end

function delete(client::K8sClient, resource::Resource, name::String)::Dict{String, Any}
    uri = ""
    if resource.group == "core"
        uri = joinpath("/api/", resource.version, resource.plural)
    else
        uri = joinpath("/apis/", resource.group, resource.version, resource.plural)
    end
    resp = request(client, "DELETE", uri)
    str = String(resp.body)
    jobj = JSON2.Parser.parse(str)
    return jobj
end

function delete(client::K8sClient, resource::Resource, name::String, namespace::String)::Dict{String, Any}
    if resource.group == "core"
        uri = joinpath("/api/", resource.version, "namespaces", namespace, resource.plural)
    else
        uri = joinpath("/apis/", resource.group, resource.version, "namespaces", namespace, resource.plural)
    end
    resp = request(client, "DELETE", uri)
    str = String(resp.body)
    jobj = JSON2.Parser.parse(str)
    return jobj
end

max_stars_repo_name: pbarker/kubemeta

max_issues_repo_name: pbarker/kubemeta

max_forks_repo_name: pbarker/kubemeta


text:
using Gadfly

set_default_plot_size(6inch, 3inch)

plot(xgroup=["A", "A", "B", "B"],
     x=["X", "Y", "X", "Y"],
     y=rand(4),
     Geom.subplot_grid(Geom.bar))

plot(x=rand(3),
     y=rand(3),
     xgroup=["A", "A", "B"],
     ygroup=["X", "Y", "Y"],
     Geom.subplot_grid(Geom.bar))


text:
using TrussMorph
using Test

fp = "/Users/yijiangh/.julia/dev/TrussMorph/test/truss_json/2D_truss.json"
truss = parse_truss_json(fp)

load_fp = "/Users/yijiangh/.julia/dev/TrussMorph/test/truss_json/2D_truss_load_case.json"
load = parse_load_json(load_fp)

@test size(truss.X,2) == 2
@test size(truss.T,2) == 2
@test size(truss.S,2) == 4
@test size(load, 2) == 4
# @test_throws MethodError parse_truss_json(fp)


text:
using AMQPClient
const VIRTUALHOST = "/"
const HOST = "127.0.0.1"

function receive()
    # 1. Create a connection to the localhost or 127.0.0.1 of virtualhost '/'
    connection(; virtualhost=VIRTUALHOST, host=HOST) do conn
        # 2. Create a channel to send messages
        channel(conn, AMQPClient.UNUSED_CHANNEL, true) do chan
            # 3. Declare a exchange
            exchange = "logs"
            exchange_declare(chan, exchange, EXCHANGE_TYPE_FANOUT)
            result, queue, _, _ = queue_declare(chan, "", durable=true)
            # 4. Bind the queue
            queue_bind(chan, queue, exchange, EXCHANGE_TYPE_FANOUT)

            println(" [*] Waiting for messages. To exit press CTRL+C")
            on_receive = (msg) -> begin
                data = String(msg.data)
                println("Received the message: $data")
                basic_ack(chan, msg.delivery_tag)
            end

            success, consumer_tag = basic_consume(chan, queue, on_receive)

            while true
                sleep(1)
            end
        end
    end
end

# Don't exit on Ctrl-C
Base.exit_on_sigint(false)
try
    receive()
catch ex
    if ex isa InterruptException
        println("Interrupted")
    else
        println("Exception: $ex")
    end
end

max_forks_repo_name: tomashq/rabbitmq-tutorials


text:
# 2bit format
# ===========
#
# Reader and writer of the 2bit file format.
#
# The 2bit file format stores multiple DNA sequences (up to 4 Gbp total) as a
# compact randomly-accessible format.
# See https://genome.ucsc.edu/FAQ/FAQformat.html#format7 for the details.
#
# This file is a part of BioJulia.
# License is MIT: https://github.com/BioJulia/Bio.jl/blob/master/LICENSE.md

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