text:
using Plots; pyplot()
using Reproduce
using FileIO
using Statistics
using ProgressMeter
using JLD2
# import ..GeneralPlotUtils
# import ..LabelUtils
# include("./plot_utils.jl")
include("./plot_utils.jl")

GPU = GeneralPlotUtils
LU = LabelUtils

data_home
data_key = :twod_grid_world_error_center_dpi
# data_key = :oned_tmaze_dmu_error

folder_name = "twod_gpi"
data_home = "../TwoDGridWorld_gpi_StepSizes"

function load_data()
    experiment_folders = [data_home]
    # folder_name = "oned_rr_dpi"
    # folder_name = "oned_control_dpi"
    ic = ItemCollection(joinpath(experiment_folders[1], "data"))
    return ic
end

function load_best(ic)
# data_key = :oned_tmaze_dpi_error

    algo_divisor_keys = ["behaviour_learner", "demon_learner", "demon_opt", "demon_update"]
    sweep_params = ["demon_alpha_init", "demon_eta", "alpha_init"]

    algo_specs_full = GPU.split_algo(ic, algo_divisor_keys)

    all_algos_ics = [search(ic,algo_spec) for algo_spec in algo_specs_full]
    @show length.(all_algos_ics)

    valid_algos_ind = findall(a -> length(a) != 0, all_algos_ics)

    algo_specs = algo_specs_full[valid_algos_ind]
    algo_ics = all_algos_ics[valid_algos_ind]
    best_per_algo_ics = []
    for (i,algo_ic) in enumerate(algo_ics)
        push!(best_per_algo_ics, GPU.get_best_final_perf(algo_ic,sweep_params, data_key, 0.1))
    end
    @show length.(best_per_algo_ics)
    return best_per_algo_ics
end

ic = load_data()
best_ics = load_best(ic)

stepsizes = [GPU.load_results(algo_ic,:autostep_stepsize; return_type = "array") for algo_ic in best_ics]
#algo x runs x steps x [SF, params]

ss_algos = []
ss_algos_SF = []
for algo_ss in stepsizes
    ss_algo = []
    ss_algo_SF = []
    for run_ss in algo_ss
        # @show size(run_ss)
        # @show size(run_ss[1])
        # @show size(run_ss[1][1])
        r = []
        SF = []
        for step in run_ss
            if size(step[1]) == (4,4)
                push!(r, step[1])
                push!(SF,step[2])
            else
                push!(r,step[2])
                push!(SF,step[1])
            end
            # @show size(r[end])
        end
        push!(ss_algo, cat(r..., dims=3))
        push!(ss_algo_SF, cat(SF..., dims=3))
        # tmp_r = cat([step[1] for step in run_ss]..., dims = 3)
        # tmp_sf = cat([step[2] for step in run_ss]...,dims=3)

        # @show size(tmp_r)
        # push!(ss_algo, tmp_r)
    end
    # @show size(cat(ss_algo..., dims = 4))
    push!(ss_algos,cat(ss_algo..., dims = 4))
    push!(ss_algos_SF,cat(ss_algo_SF..., dims = 4))

end
function extract_cumulants_ss(ss_algo)
    distractor = ss_algo[1,1,:,:]
    drifter = ss_algo[3,3,:,:]
    c1 = ss_algo[2,2,:,:]
    c2 = ss_algo[4,4,:,:]
    return distractor,c1,drifter,c2
end

ps = [plot(title="Distractor"),plot(title="Constant"),plot(title="Drifter"),plot(title="Constant")]
ps_SF = [plot(title="Distractor"),plot(title="Constant"),plot(title="Drifter"),plot(title="Constant")]

for (ind,ss_algo) in enumerate(ss_algos)
    demon_update = best_ics[ind].items[1].parsed_args["demon_update"]
    dis,c1,drif,c2 = extract_cumulants_ss(ss_algo)

    #Rewards
    plot!(ps[1],mean(dis,dims=2), label = demon_update, legend=:topleft)
    plot!(ps[2],mean(c1,dims=2), label = demon_update, legend=:topleft)
    plot!(ps[3],mean(drif,dims=2), label = demon_update, legend=:topleft)
    plot!(ps[4],mean(c2,dims=2), label = demon_update, legend=:topleft)

    #SF
    @show size(ss_algos_SF[ind])
    dis_SF = mean(ss_algos_SF[ind][1:4,:,:,:], dims = [1,2,4])
    c1_SF = mean(ss_algos_SF[ind][5:8,:,:,:], dims = [1,2,4])
    drif_SF = mean(ss_algos_SF[ind][9:12,:,:,:], dims = [1,2,4])
    c2_SF = mean(ss_algos_SF[ind][13:16,:,:,:], dims = [1,2,4])

    @show mean(dis_SF)
    @show mean(c1_SF)
    @show mean(drif_SF)
    @show mean(c2_SF)

    @show vec(drif_SF)[end]

    plot!(ps_SF[1],vec(dis_SF), label = demon_update, legend=:topleft)
    plot!(ps_SF[2], vec(c1_SF), label = demon_update, legend=:topleft)
    plot!(ps_SF[3], vec(drif_SF), label = demon_update, legend=:topleft)
    plot!(ps_SF[4], vec(c2_SF), label = demon_update, legend=:topleft)

end
display(plot(ps...))
display(plot(ps_SF...))


text:
export WeightedRandomPolicy

using Random
using StatsBase: sample, Weights

"""
    WeightedRandomPolicy(actions, weight, sums, rng)

Similar to [`RandomPolicy`](@ref), but the probability of
each action is set in advance instead of a uniform distribution.

- `actions` are all possible actions.
- `weight` can be an 1-D or 2-D array. If it's 1-D, then the `weight` applies to all states. If it's 2-D, then the state is assume to be of `Int` and for different state, the corresponding weight is selected.
"""
struct WeightedRandomPolicy{N,A,W<:AbstractArray,S,R<:AbstractRNG} <: AbstractPolicy
    actions::A
    weight::W
    sums::S
    rng::R
end

function WeightedRandomPolicy(
    weight::W;
    actions = axes(weights, 1),
    seed = nothing,
) where {W<:AbstractArray}
    rng = MersenneTwister(seed)
    N = ndims(W)

    if N == 1
        sums = sum(weight)
    elseif N == 2
        sums = vec(sum(weight, dims = 1))
    end
    WeightedRandomPolicy{ndims(W),typeof(actions),W,typeof(sums),typeof(rng)}(
        actions,
        weight,
        sums,
        rng,
    )
end

Random.seed!(p::WeightedRandomPolicy, seed) = Random.seed!(p.rng, seed)

RLBase.update!(p::WeightedRandomPolicy, experience) = nothing

(p::WeightedRandomPolicy{1})(obs, ::MinimalActionSet) =
    sample(p.rng, p.actions, Weights(p.weight, p.sums))

function (p::WeightedRandomPolicy{1})(obs, ::FullActionSet)
    legal_actions = get_legal_actions(obs)
    legal_actions_mask = get_legal_actions_mask(obs)
    masked_weight = @view p.weight[legal_actions_mask]
    legal_actions[sample(p.rng, Weights(masked_weight))]
end

function (p::WeightedRandomPolicy{2})(obs, ::MinimalActionSet)
    s = get_state(obs)
    weight = @view p.weight[:, s]
    sample(p.rng, p.actions, Weights(weight, p.sums[s]))
end

function (p::WeightedRandomPolicy{2})(obs, ::FullActionSet)
    s = get_state(obs)
    legal_actions = get_legal_actions(obs)
    legal_actions_mask = get_legal_actions_mask(obs)
    masked_weight = @view p.weight[legal_actions_mask, s]
    legal_actions[sample(p.rng, Weights(masked_weight))]
end

ext: jl


text:
using OhMyREPL
using Test

@testset "overwriting" begin
    @testset "force overwrite" begin
        @info "A"^8 overwrite_lastlog=true
        sleep(0.2)
        @info "B"^8 overwrite_lastlog=true
        sleep(0.2)
        @info "C"^8 overwrite_lastlog=true
        sleep(0.2)
        @info "D"^8 overwrite_lastlog=true
        sleep(0.2)
        @info "E"^8 overwrite_lastlog=true
        sleep(0.2)

        @info "F"^3 n=2 overwrite_lastlog=true
        sleep(0.2)
        @info "G"^3 n=3 m=3  overwrite_lastlog=true
         sleep(0.2)
        @info "H"^3  overwrite_lastlog=true

    end

    @testset "repeat overwrite" begin
        for ii in 1:10
            @info "Should Overwrite1" ii
            sleep(0.2)
        end

        for ii in 1:10
            @info "Should Overwrite2 (but not first)" ii doubleii=2ii
            sleep(0.2)
        end
    end

    @testset "don't repeat overwrite" begin
        for ii in 1:10
            @info "Should not overwrite" ii overwrite_lastlog=false
            sleep(0.2)
        end
    end
end

@testset "progress" begin

    @testset "single" begin
        @info "A" progress=0.1
        @info "A" progress=-0.1
        @info "A" progress=10
        @info "A" progress=NaN
        @info "A" progress=Inf
        @info "A" progress=-Inf

        @info "A" progress="Good"

        @info "B" progress = "50%"
        @info "B" progress = "fully done 100%"

        @info "C" completion = "5%"
    end

    @testset "loop" begin
        for ii in 1:100
            @info "loop" ii progress=ii/100 half="$(ii/2)%"
            sleep(0.02)
        end
    end

end


text:
using DelimitedFiles, Memoize;

input = readdlm(ARGS[1], ')', String, '\n');

# The orbits dictionary is a mapping of moon -> planet -> star -> etc.
# A relation in this set means "orbits" ("(").
orbits = Dict();
foreach(i -> orbits[input[i,2]] = input[i,1], 1:size(input)[1])

# Count the orbits by exploring the graph until we reach the root.
# We can improve this by memoizing each call.
@memoize function countOrbits(object)
    #println("$(object) in orbits: $(haskey(orbits, object))")
    if haskey(orbits, object)
        return 1 + countOrbits(orbits[object])
    else
        return 0
    end
end

println("Day 6 Part 1: $(reduce(+, map(countOrbits, collect(keys(orbits)))))")
println("Day 6 Part 2:\n  Reduce your input with the following command and solve at https://wjholden.com/dijkstra")
println("""  Get-Content $(ARGS[1]) | ForEach-Object { "\$_".Replace(")", " ") + " 1"; }""")
println("  The correct answer is 2 less than the computed distance.")


text:
module FloatProps

export FloatTypes, IntTypes
export floatType, intType
export mkFloat, mkInt, asFloat, asInt
export float_bits
export mantissa_bits, signbit_bits, exponent_bits
export mantissa_mask, exponent_mask, signbit_mask
export exponent_offset, min_exponent, max_exponent

make_mask{T}(nbits::T, nshift::T) =
    convert(T, (convert(T,1) << nbits - convert(T,1)) << nshift)

# Defined in Base:
# eps, nextfloat, precision, prevfloat, realmax, realmin, typemax, typemin
# Inf, Nan

typealias FloatTypes Union(Float16, Float32, Float64)
typealias IntTypes Union(Int16, Int32, Int64)

@inline floatType{T<:FloatTypes}(::Type{T}) = T
@inline intType(::Type{Float16}) = Int16
@inline intType(::Type{Float32}) = Int32
@inline intType(::Type{Float64}) = Int64

@inline mkFloat{T<:FloatTypes}(::Type{T}, x) = convert(T, x)
@inline mkInt{T<:FloatTypes}(::Type{T}, x) = convert(intType(T), x)
@inline asFloat{T<:FloatTypes}(::Type{T}, x) = reinterpret(T, x::intType(T))
@inline asInt{T<:FloatTypes}(::Type{T}, x) = reinterpret(intType(T), x::floatType(T))

@inline float_bits{T<:FloatTypes}(::Type{T}) = mkInt(T, 8*sizeof(T))

@inline mantissa_bits{T<:FloatTypes}(::Type{T}) =
    mkInt(T, precision(convert(T,0)) - 1)
@inline signbit_bits{T<:FloatTypes}(::Type{T}) =
    mkInt(T, 1)
@inline exponent_bits{T<:FloatTypes}(::Type{T}) =
    float_bits(T) - mantissa_bits(T) - signbit_bits(T)

@inline mantissa_mask{T<:FloatTypes}(::Type{T}) =
    make_mask(mantissa_bits(T), mkInt(T,0))
@inline exponent_mask{T<:FloatTypes}(::Type{T}) =
    make_mask(exponent_bits(T), mantissa_bits(T))
@inline signbit_mask{T<:FloatTypes}(::Type{T}) =
    make_mask(signbit_bits(T), exponent_bits(T) + mantissa_bits(T))

@inline exponent_offset{T<:FloatTypes}(::Type{T}) =
    mkInt(T,1) << (exponent_bits(T) - mkInt(T,1)) - mkInt(T,1)
@inline min_exponent{T<:FloatTypes}(::Type{T}) =
    mkInt(T,2) - exponent_offset(T)
@inline max_exponent{T<:FloatTypes}(::Type{T}) =
    mkInt(T,1) << exponent_bits(T) - exponent_offset(T) - mkInt(T,2)

end


text:
# SPDX-License-Identifier: X11
# 2020-08-30
# Sum of product of pairs (300pt)

const M = 1000000007

function main()
  N = parse(Int, readline())

  s = parse(Int, readuntil(stdin, " "))
  t = 0
  for i ∈ 2:N
    A = parse(Int, readuntil(stdin, " "))
    t = mod(t + A * s, M)
    s = mod(s + A, M)
  end

  println(t)
end

main()

max_stars_repo_name: jthvai/competitive-programming

max_issues_repo_name: jthvai/competitive-programming

max_forks_repo_name: jthvai/competitive-programming


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 wing_example(; prompt=true, file_name="temp_wing00")

  # PARAMETERS
  n_up = 40             # Number of sections in upper surface of wing
  n_lower = 20          # Number of sections in lower surface of wing
  z1 = 0.0              # Position of first airfoil
  z2 = 10.0             # Position of second airfoil
  c1 = 2.5              # Chord of first airfoil
  c2 = 1.0              # Chord of second airfoil
  # file_name = "temp_wing00"    # Name for vtk outpout file

  # Reads the original airfoil geometry from airfoiltools.com
  org_x1, org_y1 = readcontour(data_path*"naca6412.dat"; header_len=1)
  org_x2, org_y2 = readcontour(data_path*"naca0006.dat"; header_len=1)

  # Separate upper and lower sides to make the contour injective in x
  upper1, lower1 = splitcontour(org_x1, org_y1)
  upper2, lower2 = splitcontour(org_x2, org_y2)

  # Parameterize both sides independently
  fun_upper1 = gt.parameterize(upper1[1], upper1[2], zeros(eltype(upper1[1]), size(upper1[1])); inj_var=1)
  fun_lower1 = gt.parameterize(lower1[1], lower1[2], zeros(eltype(lower1[1]), size(lower1[1])); inj_var=1)
  fun_upper2 = gt.parameterize(upper2[1], upper2[2], zeros(eltype(upper2[1]), size(upper2[1])); inj_var=1)
  fun_lower2 = gt.parameterize(lower2[1], lower2[2], zeros(eltype(lower2[1]), size(lower2[1])); inj_var=1)

  # Upper surface sections
  aux1 = Int(floor(20/53*n_up))
  sec1 = (0.35, aux1, 3.0, true) # 35% of the line has 20 sections in ratio 3.0 around center
  sec2 = (0.65, n_up-aux1, 3.0, true) # 65% of the line has 33 sections in ratio 3.0 around center

  # New discretization for both surfaces
  upper_points1 = gt.multidiscretize(fun_upper1, 0, 1, [sec1,sec2])
  lower_points1 = gt.discretize(fun_lower1, 0, 1, n_lower, 8.0; central=true)
  upper_points2 = gt.multidiscretize(fun_upper2, 0, 1, [sec1,sec2])
  lower_points2 = gt.discretize(fun_lower2, 0, 1, n_lower, 8.0; central=true)

  # Put both surfaces back together from TE over the top and from LE over the bottom.
  reverse!(upper_points1)                           # Trailing edge over the top
  new_x1 = [point[1] for point in upper_points1]
  new_y1 = [point[2] for point in upper_points1]     # Leading edge over the bottom
  new_x1 = vcat(new_x1, [point[1] for point in lower_points1])
  new_y1 = vcat(new_y1, [point[2] for point in lower_points1])
  reverse!(upper_points2)                           # Trailing edge over the top
  new_x2 = [point[1] for point in upper_points2]
  new_y2 = [point[2] for point in upper_points2]     # Leading edge over the bottom
  new_x2 = vcat(new_x2, [point[1] for point in lower_points2])
  new_y2 = vcat(new_y2, [point[2] for point in lower_points2])

  plot_airfoil(new_x1, new_y1; style="--.k", title_str="NACA6412")
  println("Close figure and press ENTER")
  readline()
  plot_airfoil(new_x2, new_y2; style="--.k", title_str="NACA0006")
  println("Close figure and press ENTER")
  readline()

  npoints = size(new_x1)[1]         # Number of points on each airfoil contour
  # Dimensions the airfoil acording to their chord length and positions them in z
  airfoil1 = [[c1*new_x1[i], c1*new_y1[i], z1] for i in 1:npoints]
  airfoil2 = [[c2*new_x2[i], c2*new_y2[i], z2] for i in 1:npoints]

  # Discretization of the lofting between the airfoils
  nloft = 50                      # Number of sections
  rloft = 10.0                   # Expansion ratio
  central = true                  # Expands about center between edges
  sections = [(1.0, nloft, rloft, central)]   # Discretization sections

  # Dummy point data for good looking visuals
  pd1 = [i for i in 1:size(airfoil1)[1]]
  pd2 = size(airfoil1)[1] .+ [i for i in 1:size(airfoil2)[1]]

  # Generates cells in VTK Legacy format
  out = gt.lines2vtkmulticells(airfoil1, airfoil2, sections;
                                          point_data1=pd1, point_data2=pd2)
  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)

  # 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:
# Gibbs samplers for Beta NTL models of partitions, graphs
using Distributions
using StatsBase

include("../dataset.jl")
include("../crp.jl")
include("../ntl_gibbs.jl")
include("../slice.jl")
include("../evaluation.jl")
include("../gibbs_util.jl")
include("../geometric_ia.jl")
include("../poisson_ia.jl")

###########################################################################
# 1. Gibbs sampler settings
###########################################################################
dataset_name = "synthetic geometric" # set dataset; see section 2. Data for options that are automated

arrivals = "geometric" # set arrival time model
save_output = false # whether or not to save sampler output

if startswith(dataset_name, "synthetic")
  N = 1000 # set size of synthetic data (# of ends of edges)
end

n_iter = 10000 # 50000  # total number of Gibbs iterations to run
n_burn = 5000   # burn-in
n_thin = 100     # collect every `n_thin` samples

n_print = 1000 # prints updates every `n_print` iterations

# set which components to update
gibbs_psi = true            # NTL Ψ paramters
gibbs_alpha = true          # NTL alpha parameter
gibbs_arrival_times = true  # arrival times
gibbs_ia_params = true     # arrival time distribution parameters
gibbs_perm_order = true   # order of blocks in partition/vertices in graph

## SET SEED
srand(0)

############################################################################
# 2. Data
############################################################################

if startswith(dataset_name,"synthetic crp") # Synthetic data w/ CRP interarrivals
  println("Synethsizing data.")
  # set generating parameter values
  ntl_alpha = 0.8
  crp_theta = 1.0
  if endswith(dataset_name,"uncoupled")
    crp_alpha = 0.6
  else
    crp_alpha = ntl_alpha
  end

  # create intearrival distribution object and synthetic data
  interarrival_dist = CRP(crp_theta,crp_alpha)
  Z_data, PP_data, T_data = generateLabelSequence(N,ntl_alpha,interarrival_dist)

  # function check
  assert(all(PP_data .== seq2part(Z_data)))
  assert(all(T_data .== get_arrivals(Z_data)))

elseif dataset_name=="synthetic geometric" # Synthetic data w/ geometric interarrivals
  println("Synethsizing data.")
  ntl_alpha = 0.8
  geom_p = 0.5

  # create intearrival distribution object and synthetic data
  interarrival_dist = Geometric
  Z_data, PP_data, T_data = generateLabelSequence(N,ntl_alpha,interarrival_dist(geom_p))

  # function check
  assert(all(PP_data .== seq2part(Z_data)))
  assert(all(T_data .== get_arrivals(Z_data)))

elseif dataset_name=="synthetic poisson" # Synthetic data w/ Poisson interarrivals
  println("Synethsizing data.")
  ntl_alpha = 0.5
  poisson_lambda = 4.0

  # create intearrival distribution object and synthetic data
  interarrival_dist = Poisson
  Z_data, PP_data, T_data = generateLabelSequence(N,ntl_alpha,interarrival_dist(poisson_lambda))

  # function check
  assert(all(PP_data .== seq2part(Z_data)))
  assert(all(T_data .== get_arrivals(Z_data)))

elseif dataset_name=="college msg" # use College Message data
  elist = readdlm("../data/CollegeMsg.txt",Int64)
  Z_data = vec(elist[:,1:2]')
  PP_data = seq2part(Z_data)
  T_data = get_arrivals(Z_data)

elseif dataset_name=="mathoverflow"
  Z_data = vec(readdlm("../data/sx-mathoverflow-Z.txt",Int64))
  PP_data = seq2part(Z_data)
  T_data = get_arrivals(Z_data)

end

# sort partition if necessary
if gibbs_perm_order
  PP_sort = sortrows(hcat(PP_data,collect(1:size(PP_data,1))),rev=true)
  # sort ties in ascending order of original order for plotting purposes
  maxdeg = maximum(PP_sort)
  for j in 1:maxdeg
    PP_sort[PP_sort[:,1].==j,2] = sort(PP_sort[PP_sort[:,1].==j,2],rev=false)
  end
  perm_data = PP_sort[:,2]
  PP = PP_sort[:,1]

else
  PP = deepcopy(PP_data)
  perm_data = collect(1:size(PP,1))

end

K = size(PP,1)
N = sum(PP)
gibbs_alpha ? nothing : alpha_fixed = ntl_alpha

println("Finished pre-processing data.")

###########################################################################
# 3. Settings for sampling the BNTL α parameter
###########################################################################

# Gamma prior parameters
ntl_gamma_a = 1.
ntl_gamma_b = 10.

# prior is specified as a distribution on (0,Inf);
# transformation to α ∈ (-Inf,1) will be performed during sampling as appropriate
ntl_alpha_prior_dist = Gamma(ntl_gamma_a,ntl_gamma_b)
ntl_alpha_log_prior = x -> logpdf(ntl_alpha_prior_dist,x)

w_alpha = 5.0 # slice sampling "window" parameter

###########################################################################
# 4. Modeled arrival time distribution settings
###########################################################################

if startswith(arrivals,"crp")
  # CRP arrival distribution
  ia_dist = v -> CRP(v[1],v[2])

  # prior on CRP parameters
  n_ia_params = 2
  theta_gamma_a = 0.1
  theta_gamma_b = 10.
  alpha_beta_a = 1.
  alpha_beta_b = 1.
  ia_prior_params = [theta_gamma_a; # prior on theta
                     theta_gamma_b;
                     alpha_beta_a; # prior on crp_alpha
                     alpha_beta_b]
 ia_theta_prior = Gamma(theta_gamma_a,theta_gamma_b)
 ia_alpha_prior = Beta(alpha_beta_a,alpha_beta_b)
 lpdf_ia_param_prior = pp -> logpdf(ia_theta_prior,pp[1]+pp[2]) + logpdf(ia_alpha_prior,pp[2])

  # slice sampling functions/parameters
  f_lp_t = (x,a) -> logpdf(ia_theta_prior,x+a)
  f_lp_a = x -> logpdf(ia_alpha_prior,x)
  w_t = 1.0 # slice sampling w parameter for crp_theta
  w_a = 1.0 # slice sampling w parameter for crp_alpha

  # CRP-specific sampling functions
  initialize_arrival_params = v -> initialize_crp_params(Gamma(v[1],v[2]),Beta(v[3],v[4]))
  if arrivals=="crp-uncoupled"
    update_arrival_params! = (ap,PP,T,n,pripar) -> update_crp_interarrival_params!(ap,PP,T,n,f_lp_t,f_lp_a,w_t,w_a,false)
  elseif arrivals=="crp-coupled"
    update_arrival_params! = (ap,PP,T,n,pripar) -> update_crp_interarrival_params!(ap,PP,T,n,f_lp_t,f_lp_a,w_t,w_a,true)
  else
    error("Unsupported arrival distribution specification.")
  end

  gibbs_ia_params ? nothing : arrival_params_fixed = [mean(Gamma(ia_prior_params[1],ia_prior_params[2])); mean(Beta(ia_prior_params[3],ia_prior_params[4]))]

elseif arrivals=="geometric"
  # Geometric interarrival distribution
  ia_dist = p -> Geometric(p[1])
  a_beta = 1.
  b_beta = 1.
  ia_prior_params = [a_beta; b_beta]
  ia_param_prior = Beta(a_beta,b_beta)
  lpdf_ia_param_prior = pp -> logpdf(ia_param_prior,pp[1])
  n_ia_params = 1
  # set update functions
  gibbs_ia_params ? nothing : arrival_params_fixed = [mean(ia_param_prior)]
  initialize_arrival_params = v -> [(K-1)/(N-1)]
  update_arrival_params! = update_geometric_interarrival_param!

elseif arrivals=="poisson"
  # Poisson interarrival distribution
  ia_dist = lambda -> Poisson(lambda[1])
  a_gamma = 0.01
  b_gamma = 10.
  ia_prior_params = [a_gamma; b_gamma]
  ia_param_prior = Gamma(a_gamma,b_gamma)
  lpdf_ia_param_prior = pp -> logpdf(ia_param_prior,pp[1])
  n_ia_params = 1
  # set update functions
  gibbs_ia_params ? nothing : arrival_params_fixed = mean(ia_param_prior)
  initialize_arrival_params = v -> [(N-1)/(K-1)]
  update_arrival_params! = update_poisson_interarrival_param!

else
  error("Unsupported arrival distribution specification.")
end

############################################################################
# 5. Gibbs sampler
############################################################################

# PP is the vector of block sizes/vertex degrees to be used.
# If block/vertex order is being inferred, PP should be sorted in descending order;
#   otherwise, it should be in arrival-order

sampler_control = GibbsSamplerControls(
                    n_iter,n_burn,n_thin,n_print,
                    gibbs_psi,gibbs_alpha,gibbs_arrival_times,gibbs_ia_params,gibbs_perm_order,
                    arrivals,ia_dist,ia_prior_params,n_ia_params,initialize_arrival_params,update_arrival_params!,
                    ntl_alpha_log_prior,ntl_alpha_prior_dist,lpdf_ia_param_prior,
                    w_alpha,
                    0.5,[0.5] # these are place-holders that aren't used when all parameters are being inferred
                  )

spl_out = GibbsSampler(sampler_control,PP,T_data)

############################################################################
# 6. Save sample output
############################################################################

if save_output
  using JLD
  using JSON
  using DataStructures

  datetime = Dates.format(now(), "yyyy-mm-dd_HH-MM-SS")
  try
    Base.Filesystem.mkdir("./sampler_output")
  catch
    nothing
  end
  dirname = "./sampler_output/gibbs_" * datetime * "_"
  fname = "samples.jld"
  pathname = dirname * fname
  save(pathname,
      "psi_gibbs",spl_out.psi,
      "T_gibbs",spl_out.T,
      "alpha_gibbs",spl_out.alpha,
      "ia_params_gibbs",spl_out.ia_params,
      "perm_gibbs",spl_out.sigma,
      "N",N,"K",K,"t_elapsed",spl_out.t_elapsed,
      "PP",PP,"perm_data",perm_data,"T_data",T_data,"Z_data",Z_data)

  params = OrderedDict(
  "dataset_name" => dataset_name, "arrivals" => arrivals,
  "n_iter" => n_iter, "n_burn" => n_burn, "n_thin" => n_thin,
  "gibbs_psi" => gibbs_psi, "gibbs_alpha" => gibbs_alpha, "gibbs_arrival_times" => gibbs_arrival_times, "gibbs_ia_params" => gibbs_ia_params, "gibbs_perm_order" => gibbs_perm_order,
  "ntl_gamma_a" => ntl_gamma_a, "ntl_gamma_b" => ntl_gamma_b, "w_alpha" => w_alpha
  )
  if dataset_name == "synthetic crp"
      params["N"] = N
      params["ntl_alpha"] = ntl_alpha
      params["crp_theta"] = crp_theta
      params["crp_alpha"] = crp_alpha
  elseif dataset_name == "synthetic geometric"
      params["N"] = N
      params["ntl_alpha"] = ntl_alpha
      params["geom_p"] = geom_p
  end
  if arrivals == "crp"
      params["theta_gamma_a"] = theta_gamma_a
      params["theta_gamma_b"] = theta_gamma_b
      params["alpha_beta_a"] = alpha_beta_a
      params["alpha_beta_b"] = alpha_beta_b
      params["w_t"] = w_t
      params["w_a"] = w_a
  elseif arrivals == "geometric"
      params["a_beta"] = a_beta
      params["b_beta"] = b_beta
  end

  open(dirname * "params.json", "w") do f
      write(f, JSON.json(params))
  end
end


text:
using YaoArrayRegister
using Test

@testset "test utils" begin
    include("utils.jl")
end

@testset "test ArrayReg" begin
    include("register.jl")
    include("operations.jl")
    include("focus.jl")
end

@testset "test instructions" begin
    include("instruct.jl")
end

@testset "test measure" begin
    include("measure.jl")
end


text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using Base.BinaryPlatforms
include("../common.jl")

gap_version = v"400.1191.001"
gap_lib_version = v"400.1191.000"
name = "JuliaInterface"
upstream_version = v"0.7.3" # when you increment this, reset offset to v"0.0.0"
offset = v"0.0.0" # increment this when rebuilding with unchanged upstream_version, e.g. gap_version changes
version = offset_version(upstream_version, offset)

julia_versions = [v"1.6.0", v"1.7.0", v"1.8.0"]

# Collection of sources required to build libsingular-julia
sources = [
    GitSource("https://github.com/oscar-system/GAP.jl", "baa0589573a9b56a01d850c1c0b5a381fbec07bb"),
]

# Bash recipe for building across all platforms
script = raw"""
cd GAP.jl/pkg/JuliaInterface
./configure --with-gaproot=${prefix}/share/gap
make -j${nproc} CFLAGS="-I${includedir}"

# copy the loadable module
mkdir -p ${prefix}/lib/gap
cp bin/*/*.so ${prefix}/lib/gap/

# copy the sources, too, so that we can later compare them
cp -r src ${prefix}/

install_license ../../LICENSE
"""

name = gap_pkg_name(name)
platforms, dependencies = setup_gap_package(gap_version, gap_lib_version)

# expand julia platforms
julia_platforms = []
for p in platforms
    for jv in julia_versions
        if jv == v"1.6.0" && Sys.isapple(p) && arch(p) == "aarch64"
            continue
        end
        p = deepcopy(p)
        BinaryPlatforms.add_tag!(p.tags, "julia_version", string(jv))
        push!(julia_platforms, p)
    end
end

push!(dependencies, BuildDependency("libjulia_jll"))

# The products that we will ensure are always built
products = [
    FileProduct("lib/gap/JuliaInterface.so", :JuliaInterface),
]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, julia_platforms, products, dependencies;
               julia_compat="1.6", preferred_gcc_version=v"7")

max_stars_repo_name: sharanry/Yggdrasil

max_issues_repo_name: sharanry/Yggdrasil

max_forks_repo_name: sharanry/Yggdrasil


text:
# ======================== mendeleev(Z) ========================================

@doc raw"""
    mendeleev(Z::Int)

The properties `name` and `symbol` of the *element* with *atomic number* `Z`.
#### Example:
```
mendeleev(11)
 ("Sodium", "Na")
```
"""
function mendeleev(Z::Int)

    0 < Z < 45 || return error("jwError: element not inplemeted")

    element = Dict(1 => ("Hydrogen", "H"), 2 => ("Helium", "He"),
    3 => ("Lithium", "Li"), 4 => ("Beryllium", "Be"), 5 => ("Boron", "B"), 6 => ("Carbon", "C"),
    7 => ("Oxygen", "O"), 8 => ("Nitrogen", "N"), 9 => ("Fluorine", "F"), 10 => ("Neon", "Ne"),
    11 => ("Sodium", "Na"), 12 => ("Magnesium", "Mg"), 13 => ("Aluminium", "Al"), 14 => ("Silicon", "Si"),
    15 => ("Phosphorus", "P"), 16 => ("Sulphor", "S"), 17 => ("Chlorine", "Cl"), 18 => ("Argon", "Ar"),
    19 => ("Potassium", "K"), 20 => ("Calcium", "Ca"), 21 => ("Scandium", "Sc"), 22 => ("Titanium", "Ti"),
    23 => ("Vanadium", "Va"), 24 => ("Chromium", "Cr"), 25 => ("Manganese", "Cl"),26 => ("Iron", "Fe"),
    27 => ("Cobalt", "Co"), 28 => ("Nickel", "Ni"), 29 => ("Copper", "Cu"), 30 => ("Zinc", "Zn"),
    31 => ("Gallium", "Ga"), 32 => ("Germanium", "Ge"), 33 => ("Arsenic", "As"), 34 => ("Selenium", "Se"),
    35 => ("Bromine", "Br"), 36 => ("Krypton", "Mg"),
    37 => ("Rubidium", "Rb"), 38 => ("Strontium", "Sr"), 39 => ("Yttrium", "Y"), 40 => ("Zirconium", "Zr"),
    41 => ("Niobium", "Nb"), 42 => ("Molybdenium", "Mo"), 43 => ("Technetium", "Tc"), 44 => ("Ruthenium", "Ru"))

    return (name,symbol) = get(element,Z,nothing)

end

# ======================== Atom(name, symbol, Z, I, Q, M, I, gI) ===============

@doc raw"""
    Atom(name::String, symbol::String, Z::Int, Q::Int, M::Float64, I::Real, gI::Float64)

Type to specify the *atomic species* with fields `name::String` (name of element),
`symbol::String` (symbol of element), `Z::Int` (atomic number),
`Q::Int` (ionic charge in a.u.), `M::Float6` (nuclear mass in amu),
`I::Rational{Int}` (nuclear spin as integer or rational number), `gI::Float64` (nuclear g-factor)

Note: the type `Atom` is best created by the function `createAtom`.
#### Examples:
```
Hydrogen = Atom("Hydrogen", "¹H", 1, 0, 1.00782503223, 1//2, 5.585694713)

name = Hydrogen.name
symbol = Hydrogen.symbol
Z = Hydrogen.Z
Q = Hydrogen.Q
M = Hydrogen.M
I = Hydrogen.I
gI = Hydrogen.gI
println("name = $(name), symbol = $(symbol), Z = $Z, Q = $Q, M = $M, I = $I, gI = $gI")
 name = Hydrogen, symbol = ¹H, Z = 1, Q = 0, M = 1.00782503223, I = 1//2, gI = 5.585694713
```
"""
struct Atom              # atomic properties
    name::String         # name of element
    symbol::String       # atomic symbol
    Z::Int               # Z: atomic number
    Q::Int               # Q: ionic charge (a.u.)
    M::Float64           # M: nuclear mass (amu)
    I::Real              # I: nuclear spin as integer or rational number
    gI::Float64          # gI nucear g-number
end

# ======================== createAtom(Z, Q, M, I, gI) ===========

@doc raw"""
    createAtom(Z::Int, Q::Int, M::Float64, I::Real, gI::Float64)

Create Atom Type with fields `name::String` (name of element), `symbol::String` (symbol of element),
`Z::Int` (atomic number), `Q::Int` (ionic charge in a.u.), `M::Float6` (nuclear mass in amu),
`I::Real` (nuclear spin as integer or rational number), `gI::Float64` (nuclear g-factor).
#### Examples:
```
createAtom(1,0,1.00782503223,1//2,5.585694713)
 Atom created: Hydrogen, symbol = ¹H, Z = 1, Q = 0, M = 1.00782503223, I = 1//2, gI = 5.585694713
 Atom("Hydrogen", "¹H", 1, 0, 1.00782503223, 1//2, 5.585694713)

createAtom(2,1,4.00260325413,1//2,0.0)
 Atom created: Helium ion, symbol = ⁴Heᐩ, Z = 2, Q = 1, M = 4.00260325413, I = 1//2, gI = 0.0
 Atom("Helium ion", "⁴Heᐩ", 2, 1, 4.00260325413, 1//2, 0.0)
```
"""
function createAtom(Z::Int, Q::Int, M::Float64, I::Real, gI::Float64)

    S = typeof(I) ∈ [Float16,Float32,Float64] ? rationalize(I) : I

    strQ = abs(Q) > 1 ? sup(abs(Q)) : ""
    strQ = Q > 0 ? (strQ * 'ᐩ') : Q < 0 ? (strQ * 'ᐨ') : ""

    (name,symbol) = mendeleev(Z)

    name = Q ≠ 0 ? (name * " ion") : name
    symbol = sup(Int(round(M))) * symbol * strQ

    println("Atom created: $(name), symbol = $(symbol), Z = $Z, Q = $Q, M = $M, I = $I, gI = $gI")

    return Atom(name, symbol, Z, Q, M, I, gI)

end

# ======================== Term(name, n, ℓ, S, L, J) ===========

@doc raw"""
    Term(name::String, n::Int, ℓ::Int, S::Real, L::Int, J::Real)

Type to specify the *fine-structure Term* in *Russell-Saunders notation* with fields
`name::String`, `n::Int` (principal quantum number),
`ℓ::Int` (orbital angular momentum valence electron),
`S::Real` (total electron spin as integer or rational number),
`L::Int` (total orbital angular momentum),
`J::Real` (total electronic angular momentum as integer or rational number).

Note: the type `Term` is best created by the function `createTerm`.
#### Examples:
```
Term_H1I = Term("1s ²S₁⸝₂", 1, 0, 1//2, 0, 1//2)
 Term("1s ²S₁⸝₂", 1, 0, 1//2, 0, 1//2)

name = Term_H1I.name
n = Term_H1I.n
ℓ = Term_H1I.ℓ
S = Term_H1I.S
L = Term_H1I.L
J = Term_H1I.J
println("name = $(name, n = $n, ℓ = $ℓ, S = $S, L = $L, J = $J")
 name = "1s ²S₁⸝₂", n = 1, ℓ = 0, S = 1//2, L = 0, J = 1//2
```
"""
struct Term
    name::String         # LS term notation
    n::Int               # principal quantum number
    ℓ::Int               # orbital angular momentum valence electron
    S::Real              # total electron spin as integer or rational number
    L::Int               # total orbital angular momentum
    J::Real              # total electronic angular momentum as integer or rational number
end

# ======================== createTerm(n, ℓ, S, L, J) ===========

@doc raw"""
    createTerm(n::Int, ℓ::Int, S::Real, L::Int, J::Real)

Specify Term Type in the *Term notatation* with fields
`n::Int` (principal quantum number), `ℓ::Int` (orbital angular momentum valence electron),
`S::Rational` (total electron spin), `L::Int` (total orbital angular momentum),
`J::Rational` (total electronic angular momentum).
#### Examples:
```
term_H1I = createTerm(1, 0, 1//2, 0, 1//2)
 Term created: 1s ²S₁⸝₂, n = 1, ℓ = 0, S = 1//2, L = 0, J = 1//2
 Term("1s ²S₁⸝₂", 1, 0, 1//2, 0, 1//2)
```
"""
function createTerm(n::Int, ℓ::Int, S::Real, L::Int, J::Real)

    S = typeof(S) ∈ [Float16,Float32,Float64] ? rationalize(S) : S
    J = typeof(J) ∈ [Float16,Float32,Float64] ? rationalize(J) : J

    strL = ['s','p','d','f','g','h','i','k','l','m','n','o','q','r','t','u']

    name = string(n) * strL[ℓ + 1] * ' ' * sup(Int(2S + 1)) * uppercase(strL[L + 1]) * sub(J)

    ℓ < n || return error("jwError: ℓ < n rule not satisfied")
    abs(L-S) ≤ J ≤ (L+S)  || return error("jwError: Δ(LSJ) condition not satisfied")

    println("Term created: $(name); n = $n, ℓ = $ℓ, S = $S, L = $L, J = $J")

    return Term(name, n, ℓ, S, L, J)

end

# ======================== Grid(N, h, r0, r, r′) ===============

@doc raw"""
    Grid(N::Int, h::Float64, r0::Float64, r::Vector{Float64}, r′::Vector{Float64})

Type to specify the `Grid` on which the atomic wavefunction is defined, with fields `N::Int` (number of grid points),
`h::Float64` (step size on uniform grid), `r0::Float64` (scale factor for physical grid in a.u.),
`r::Vector{Float64}` (tabulated grid function), r′::Vector{Float64} (tabulated derivative of grid function).

Note: the type `Grid` is created by the function `createGrid` which serves to tabulate the grid functions.
"""
struct Grid
    N::Int               # number of grid points
    h::Float64           # stepsize on uniform grid
    r0::Float64          # scale factor for physical grid in a.u.
    r::Vector{Float64}   # tabulated grid function
    r′::Vector{Float64}  # tabulated derivative of grid function
end

# ======================== gridfunction(n, h, r0; pmax=6, deriv=0)  ===============

@doc raw"""
    gridfunction(n::Int, h::Float64, r0::Float64; pmax=6, deriv=0)

Transformation from the uniform grid (used for finite-difference calculus) and the nonuniform physical grid,
```math
    f[t,h,r0]=r_0[t+\frac{1}{2!}t^2+\frac{1}{3!}t^3+\frac{1}{4!}t^4+\frac{1}{5!}t^5+\frac{1}{6!}t^6]
```
```math
    f^{′}[t,h,r0]=f[t,h,r_0;deriv=1]=r_0h[1+t+\frac{1}{2!}t^2+\frac{1}{3!}t^3+\frac{1}{4!}t^4+\frac{1}{5!}t^5],
```
where ``t=(n-1)h``.
"""
function gridfunction(n::Int, h::Float64, r0::Float64; pmax=6, deriv=0)
# =================================================================
# my grid function - exponential at small r and algebraic a large r
# gridfunction(n::Int; r0=grid.r0) = r0*(exp((n-1) * h)-1.0) # gridfunction used by Walter Johnson type
# =================================================================

    deriv < pmax || return 0.0

    t = (n-1) * h
    o = deriv > 0 ? 1.0 : 0.0
    v = 1.0

    for p=1:(pmax-deriv)
        v = v * t / p
        o = o + v
    end

    return deriv > 0 ? o * r0 * h : o * r0

end

# ======================== createGrid(N; h=0.01, r0=0.001)   ===============

@doc raw"""
    createGrid(N::Int; h=0.01, r0=0.001)

Tabulate the gridfunction for an array of `N` points with uniform grid spacing `h` and physical scale factor `r0`.
#### Example:
```
grid = createGrid(3; h=0.01, r0=0.001)
grid.r
3-element Vector{Float64}:
 0.0
 1.0050167084168057e-5
 2.020134002675555e-5
```
"""
function createGrid(N::Int; h=0.01, r0=0.001)

    r = [gridfunction(n, h, r0) for n=1:N]
    r′= [gridfunction(n, h ,r0; deriv=1) for n=1:N]

    return Grid(N, h, r0, r, r′)

end

# ============================== matG(n, Etot, atom, grid, scr) ================

@doc raw"""
    matG(n::Int, E::Float64, atom::Atom, grid::Grid, scr::Vector{Float64})

Coupling matrix for a set of two coupled differential equations,
```math
    \frac{dy}{dn}[n]=G[n]\thinspace y[n]\equiv f[n],
```
where ``y=\left(\begin{array}{cc} P\\ Q  \end{array}\right)`` and
``G=\left(\begin{array}{cc} 0 & b\\ c & 0 \end{array}\right),``
"""
function matG(n::Int, E::Float64, atom::Atom, grid::Grid, scr::Vector{Float64})
# ==============================================================================
# matG - coupling matrix - Johnson (2.54)
# ==============================================================================
    G = [0.0 1.0; 0.0 0.0]

    Z = atom.Z
    ℓ = term.ℓ
    r = grid.r[n+1]
    r′= grid.r′[n+1]
    s = scr[n+1]

    a = 0.0
    b = 1.0
    c = fGc(n, r, E, Z, ℓ, s)
    d = 0.0

    G[1,1] = a
    G[1,2] = b * r′
    G[2,1] = c * r′
    G[2,2] = d

    return G

end


text:
include("../examples/lorenz63.jl")
include("../src/clvs.jl")
#include("../src/lss.jl")
include("../src/adjoint_lss.jl")
using Test
using DiffEqSensitivity, OrdinaryDiffEq, Zygote
function test_dlorenz63()
    s = [10., 28., 8/3]
    m = 10
    u0 = rand(3,m)
    u_trj = lorenz63(u0, s, 1)[2,:,:]
    eps = 1.e-4
    du_x = (lorenz63(u_trj .+ eps.*[1; 0.; 0.],
    					s, 1)[2,:,:] -
    			lorenz63(u_trj .- eps.*[1; 0.; 0.],
    					s, 1)[2,:,:])./(2*eps)
    du_y = (lorenz63(u_trj .+ eps.*[0; 1.; 0.],
    					s, 1)[2,:,:] -
    			lorenz63(u_trj .- eps.*[0; 1.; 0.],
    					s, 1)[2,:,:])./(2*eps)
    du_z = (lorenz63(u_trj .+ eps.*[0; 0.; 1.],
    					s, 1)[2,:,:] -
    			lorenz63(u_trj .- eps.*[0; 0.; 1.],
    					s, 1)[2,:,:])./(2*eps)
    du_fd = reshape(collect([du_x; du_y; du_z]),
    					3, 3, m)
    du = dlorenz63(u_trj', s)
    @test all(isapprox.(du_fd - du, 0.,atol=1.e-8)) == 1
end
function test_les()
    s = [10., 28., 8/3]
    m = 1
    n = 20000
    u0 = rand(3,m)
    u_trj = lorenz63(u0, s, n)[:,:,1]
    du_trj = dlorenz63(u_trj, s)
    les, clv_trj = clvs(du_trj,3)
    println(les/0.005)
    @test all(isapprox.(les/0.005, [0.91, 0.,
    	-14.572], rtol=1.e-1)) == 1
end
function test_lss()
    s = [10., 28., 8/3]
    m = 1
    n = 3000
    n_runup = 5000
    u0 = rand(3,m)
    u_init = lorenz63(u0, s, n)[end,:,:]
    d = 3
    d_u = 2
    dJ = zeros(1,d,n+1)
    dJ[1,3,:] .= 1.
    n_samples = 100
    dJds = zeros(1,n_samples)
    vsh = zeros(d,n+1,n_samples)
    for i=1:n_samples
    	println("Starting LSS, sample ", i)
    	u_trj = lorenz63(u_init, s, n)[:,:,1]
    	du_trj = dlorenz63(u_trj, s)
        X = perturbation(u_trj,s) #ith col in T_{u_{i+1}} M
    	f = vectorField(u_trj,s)
    	#f = zeros(d,n+1)
    	y, xi = lss(du_trj, X, f, s, d_u)
		dJds[:,i] = compute_sens(y, xi, dJ, f)
    	@show dJds[i]
    	vsh[:,:,i] = y
    	u_init .= reshape(u_trj[end,:],3,1)
    end
    @test isapprox((sum(dJds)/n_samples),1.0,rtol=0.1)
    return vsh, dJds
end
function test_condition_number()
    s = [10., 28., 8/3]
    m = 1
    n_runup = 5000
    u0 = rand(3,m)
    u_init = lorenz63(u0, s, n_runup)[end,:,:]
    d = 3
    d_u = 2
    n_samples = 30
    dJds = zeros(n_samples)
    condnum = zeros(n_samples)
    n_arr = LinRange(500, 5000, n_samples)
    for i=1:n_samples
    	println("Starting LSS, sample ", i)
    	n = floor(Int64,n_arr[i])
    	u_trj = lorenz63(u_init, s, n)[:,:,1]
    	du_trj = dlorenz63(u_trj, s)
        X = perturbation(u_trj,s) #ith col in T_{u_{i+1}} M
    	f = zeros(d,n+1)
    	J = u_trj[:,3]
    	dJ = zeros(d,n+1)
    	dJ[3,:] .= 1.

    	y, dJds[i], c = lss(u_trj, du_trj, X, f, J,
    					 dJ, s, d_u)
    	condnum[i] = c
    	println(dJds[i])
    	u_init .= reshape(u_trj[end,:],3,1)
    end
    @test isapprox((sum(dJds)/n_samples),1.0,rtol=0.1)
    return dJds, condnum
end
function test_ensemble_sensitivity()
#let
    s = [10., 28., 8/3]
	du0, ds = zeros(3), zeros(3)
	n_samples = 1000
	#prob = ODEProblem(lorenz63_ad, u0, (0.,1.), s)
    for i =1:n_samples
		u0 = rand(3)
		println("sample number, ", i)
		#sol = solve(prob,Tsit5(),saveat=0.005,sensealg=QuadratureAdjoint())
    	du01, ds1 = Zygote.gradient(obj_fun, u0, s)
		@show du01, ds1
		du0 .+= du01/n_samples
		ds .+= ds1/n_samples
	end
	@test isapprox(ds[2], 1., rtol=0.1)
	return du0, ds
end
function test_tangent_ad()
	s = [10., 28., 8/3]
	du0, ds = zeros(3), zeros(3)
	n_samples = 1000
	#prob = ODEProblem(lorenz63_ad, u0, (0.,1.), s)
    for i =1:n_samples
		u0 = rand(3)
		println("sample number, ", i)
		#sol = solve(prob,Tsit5(),saveat=0.005,sensealg=QuadratureAdjoint())
    	du01, ds1 = Zygote.gradient(obj_fun, u0, s)
		@show du01, ds1
		du0 .+= du01/n_samples
		ds .+= ds1/n_samples
	end
	return du0, ds
end
#Tangent NILSS
function obj_fun_Q(v0,u0)
	sol=	solve(u0 + eps*v0, s, (0.,1.))
	return sol[:,end]
end
# To get Q, differentiate obj_fun_Q wrt v0
function obj_fun_v(u0, v0)
	sol=	solve(u0 + eps*v0, s + eps*ds0, (0.,1.))
	return sol[:,end]
end
# To get v, differentiate obj_fun_v wrt ds0

#Adjoint NILSS
function obj_fun_vT(w0,u0)
	sol = solve(u0, s, (0.,1.))
	sol[:, end] .+= eps*w0
	J = sum(sol[3,:])/n
	return J
end
#Derivative of obj_fun_vT wrt w0 to get vT
function obj_fun_QT(w0,u0)
	sol = solve(u0, s, (0.,1.))
	sol[:, end] .+= eps*w0
	return sol[:,end]
end
#Derivative of obj_fun_QT wrt w0 to get QT
function test_adjoint_lss()
    s = [10., 28., 8/3]
    m = 1
    n = 3000
    n_runup = 5000
    u0 = rand(3,m)
    u_init = lorenz63(u0, s, n)[end,:,:]
    d = 3
    d_u = 2
    dJ = zeros(d,n+1)
    dJ[3,:] .= 1.0/(n+1)
    n_samples = 100
    dJds = zeros(1,n_samples)
    vsh = zeros(d,n+1,n_samples)
    for i=1:n_samples
    	println("Starting LSS, sample ", i)
    	u_trj = lorenz63(u_init, s, n)[:,:,1]
		u_trj = reverse(u_trj, dims=1)
		duT_trj = permutedims(dlorenz63(u_trj, s),[2,1,3])
        X = perturbation(u_trj,s) #ith col in T_{u_{i+1}} M
    	g = vectorField(u_trj,s)
		u_next = lorenz63(reshape(u_trj[1,:],d,:),s,1)[end,:,1]
		g = [vectorField(reshape(u_next,:,d),s) g[:,2:n+1]]
		X = reshape(X, 1, d, n+1)
    	f = zeros(d,n+1)
    	J = u_trj[:,3]
    	y, xi = lss(duT_trj, dJ, f, s, 2, g)
    	vsh[:,:,i] = y
    	u_init .= reshape(u_trj[1,:],3,1)
		dJds[:,i] = compute_sens(y, xi, X, f)
		@show dJds[:,i]
    end
    #@test isapprox((sum(dJds)/n_samples),1.0,rtol=0.1)
    return vsh, dJds
end

max_stars_repo_name: nishaChandramoorthy/energies

max_issues_repo_name: nishaChandramoorthy/energies

max_forks_repo_name: nishaChandramoorthy/energies


text:
using Test

@testset "LazIO" begin
    include("testio.jl")
end


text:
using Graphics: @mustimplement

export GaussMeter, getGaussMeters, getGaussMeter,  getXValue, getYValue, getZValue, getXYZValues, getTemperature, getFrequency

abstract type GaussMeter <: Device end

Base.close(gauss::GaussMeter) = nothing

@mustimplement getXValue(gauss::GaussMeter)
@mustimplement getYValue(gauss::GaussMeter)
@mustimplement getZValue(gauss::GaussMeter)
@mustimplement getTemperature(gauss::GaussMeter)
@mustimplement getFrequency(gauss::GaussMeter)
@mustimplement calculateFieldError(gauss::GaussMeter, magneticField::Vector{<:Unitful.BField})

getGaussMeters(scanner::MPIScanner) = getDevices(scanner, GaussMeter)
function getGaussMeter(scanner::MPIScanner)
  gaussMeters = getGaussMeters(scanner)
  if length(gaussMeters) > 1
    error("The scanner has more than one gaussmeter device. Therefore, a single gaussmeter cannot be retrieved unambiguously.")
  else
    return gaussMeters[1]
  end
end

"""
Returns x,y, and z values and applies a coordinate transformation
"""
function getXYZValues(gauss::GaussMeter)
  values = [getXValue(gauss), getYValue(gauss), getZValue(gauss)]
  return gauss.params.coordinateTransformation*values
end

include("DummyGaussMeter.jl")
include("SimulatedGaussMeter.jl")
#include("LakeShore.jl")
include("LakeShoreF71.jl")


text:
#**************************************************************************************
# Calculations.jl
# =============== part of the GeoEfficiency.jl package.
#
# this file collects all the functions which are responsible for
# the calculations of the Geometrical Efficiency.
#
#**************************************************************************************

#------------------ consts - globals - imports -------------------

# parameters [relative and absolute precession - interation method ] are assigned in `Config.jl`.
# include("Config.jl")

#----------------------- GeoEff_Pnt -------------------------------

"""# UnExported

	GeoEff_Pnt(detector::CylDetector, aPnt::Point)::Float64

return the `geometrical efficiency` for the point source `aPnt` located on front
of the cylindrical detector `detector` face.

## Throw
*  an `NotImplementedError` if the point is out of the cylindrical detector `detector` face.
*  an `InValidGeometry` if the point location is invalide.

!!! note
    this is the base function that all other functions call directly or indirectly
    to calculate `geometrical efficiency` of the cylindrical-ish detector family.

"""
function GeoEff_Pnt(detector::CylDetector, aPnt::Point)::Float64
   	#aPnt.Rho > detector.CryRadius 	&& 	@notImplementedError("Point off-axis, out of the detector face")
   	detector.CryRadius > aPnt.Rho 	&& 	aPnt.Height < 0.0 	&&	@inValidGeometry("The point source location can not be inside the detector")

   	function MaxPhi(theta::Float64)::Float64
      		side = aPnt.Height * sin(theta)
      		return clamp((aPnt.Rho^2 + side^2 - detector.CryRadius^2 ) / side / aPnt.Rho / 2.0, -1.0, 1.0) |> acos
   	end # function

   	func(theta::Float64)::Float64 = MaxPhi(theta) * sin(theta)

   	if 0.0 == aPnt.Rho				# axial Point
      		strt = 0.0
      		fine = atan(detector.CryRadius, aPnt.Height)
      		return integrate(sin, strt, fine, rtol = relativeError, atol = absoluteError)[1]

   	else							# non-axial Point
      		strt = 0.0
      		transition = atan(detector.CryRadius - aPnt.Rho, aPnt.Height)
      		fine = atan(detector.CryRadius + aPnt.Rho, aPnt.Height)
      		if transition >= 0.0
	       		return integrate(sin, strt, transition, rtol = relativeError, atol = absoluteError)[1] +
                      			integrate(func, transition, fine, rtol = relativeError, atol = absoluteError)[1] / pi

      		else
				#@info("This case is not implemented yet", _file=nothing)
				@notImplementedError("Point off-axis, out of the detector face")
			# TBD: (Top + Side) efficiencies
      		end #if transition >= 0.0

   	end #if axial Point
end #function

#------------------------ GeoEff_Disk ----------------------------------

"""# UnExported

	GeoEff_Disk(detector::CylDetector, SurfacePnt::Point, SrcRadius::Real)::Float64

return the `geometrical efficiency` for a `disk` source. The `disk` center is the `SurfacePnt` and
its radius is `SrcRadius` on front of the cylindrical detector `detector` face.

produce a warning if the disk is out of the cylindrical detector face.

"""
function GeoEff_Disk(detector::CylDetector, SurfacePnt::Point, SrcRadius::Real)::Float64
   	detector.CryRadius > SurfacePnt.Rho + SrcRadius || @error(
	"off the detector face sources is not supported yet SrcRadius = $(SrcRadius), CryRadius = $(detector.CryRadius ), Rho = $(SurfacePnt.Rho)", _file=nothing)

   	integrand(xRho) = xRho * GeoEff_Pnt(detector, Point(SurfacePnt, xRho))
   	return  integrate(integrand, 0.0, SrcRadius, rtol = relativeError, atol = absoluteError)[1] / SrcRadius^2

end #function

#-------------------------- geoEff -----------------------------------

"""

	geoEff(detector::CylDetector, aSurfacePnt::Point, SrcRadius::Real = 0.0, SrcLength::Real = 0.0)::Float64

**please refer to [`geoEff(::Detector, ::Point, ::Real, ::Real)`](@ref geoEff) for more information.**

!!! warning
    `aSurfacePnt` : point `height` is considered to be measured from the detector surface.

"""
function geoEff(detector::CylDetector, aSurfacePnt::Point, SrcRadius::Real = 0.0, SrcLength::Real = 0.0)::Float64
   	detector.CryRadius < SrcRadius + aSurfacePnt.Rho   &&	@error(
		"Source Radius: Expected less than 'detector Radius=$(detector.CryRadius)', get $SrcRadius.", _file=nothing)

   	pnt::Point = deepcopy(aSurfacePnt)

   	if 0.0 == SrcRadius                         #Point source

      		detector.CryRadius > pnt.Rho || @error(
			"geoEffPoint off-axis: Expected less than 'detector Radius=$(detector.CryRadius)', get $(pnt.Rho).", _file=nothing)
        return GeoEff_Pnt(detector, pnt) / 2.0

   	elseif 0.0 == SrcLength						#Disk source

        return GeoEff_Disk(detector, pnt, SrcRadius)

   	else										# Cylindrical source

        integrand(xH::Float64) = GeoEff_Disk(detector, Point(xH, pnt.Rho), SrcRadius)
      		return integrate(integrand, aSurfacePnt.Height, aSurfacePnt.Height + SrcLength,
						rtol = relativeError, atol = absoluteError)[1] / SrcLength

   	end #if
end #function

"""

	geoEff(detector::BoreDetector, aCenterPnt::Point, SrcRadius::Real = 0.0, SrcLength::Real = 0.0)::Float64

**please refer to [`geoEff(::Detector, ::Point, ::Real, ::Real)`](@ref geoEff) for more information.**

!!! warning
    `aCenterPNT` : point `height` is consider to be measured from the detector middle, +ve value are above the detector center while -ve are below.

"""
function geoEff(detector::BoreDetector, aCenterPnt::Point, SrcRadius::Real = 0.0, SrcLength::Real = 0.0)::Float64

   	HeightWup::Float64	  = aCenterPnt.Height - detector.CryLength / 2.0
   	HeightWdown::Float64  = aCenterPnt.Height + detector.CryLength / 2.0
   	if HeightWdown < 0.0
      		if HeightWup + SrcLength < 0.0 		#invert the source.
         			return geoEff(detector, Point(aCenterPnt.Height - detector.CryLength, aCenterPnt.Rho), SrcRadius, SrcLength)

      		else # the source span the detector and emerges from both sides, split the source into two sources.
			#res = (1 - 2 * geoEff(detin, Point(0.0), SrcRadius, SrcLength))* detector.CryLength /SrcLength
         			return geoEff(detector, Point(0.0), SrcRadius, -aCenterPnt.Height) * (-aCenterPnt.Height / SrcLength) +
			       geoEff(detector, Point(0.0), SrcRadius, SrcLength + aCenterPnt.Height) * (1.0 + aCenterPnt.Height / SrcLength)

      		end
   	end

   	pntWup::Point = deepcopy(aCenterPnt);
   	aCenterPnt = Point(abs(HeightWup), aCenterPnt);  #0.0 == SrcRadius && Point(pntWup, 0.0)

   	pntWdown::Point = deepcopy(aCenterPnt);
   	pntWdown = Point(abs(HeightWdown), pntWdown); #0.0 == SrcRadius && Point(pntWdown, 0.0)

   	detin::CylDetector = CylDetector(detector.HoleRadius)
   	detout::CylDetector = CylDetector(detector.CryRadius)

   	if HeightWup >= 0.0						# the source as a whole out of detector
      		res = geoEff(detout, pntWup, SrcRadius, SrcLength) - geoEff(detin, pntWdown, SrcRadius, SrcLength)

   	elseif HeightWup + SrcLength < 0.0 		# the source as a whole in the detector
      		res = 1 - geoEff(detin, Point(abs(HeightWup + SrcLength), pntWup), SrcRadius, SrcLength)
      		res -= geoEff(detin, pntWdown, SrcRadius, SrcLength)

   	else # elseif SrcLength > 0.0
      		res = (1.0 - geoEff(detin, Point(0.0), SrcRadius, -HeightWup)) * -HeightWup / SrcLength
      		res += geoEff(detout, Point(0.0), SrcRadius, HeightWup + SrcLength) * (1.0 + HeightWup / SrcLength)

	#=else
		return 1.0 - geoEff(detin, Point(-Height, pnt), SrcRadius)[1]
	else
		res = 1 - integrate(xH -> GeoEff_Disk(detin, Point(xH, pnt), SrcRadius), 0.0, -pnt.Height, rtol=relativeError, atol=absoluteError)[1]
		res = res + integrate(xH -> GeoEff_Disk(detout, Point(xH, pntWup), SrcRadius), 0.0, pntWup.Height , rtol=relativeError, atol=absoluteError)[1]
			=#
   	end #if

   	return res
end #function

"""

	geoEff(detector::WellDetector, aWellPnt::Point, SrcRadius::Real = 0.0, SrcLength::Real = 0.0)::Float64

**please refer to [`geoEff(::Detector, ::Point, ::Real, ::Real)`](@ref geoEff) for more information.**

!!! warning
    `aWellPNT` : point `height` is considered to be measured from the detector hole surface.

"""
function geoEff(detector::WellDetector, aWellPnt::Point, SrcRadius::Real = 0.0, SrcLength::Real = 0.0)::Float64

   	pnt::Point = deepcopy(aWellPnt)
   	Height = pnt.Height - detector.HoleDepth

   	detin::CylDetector  = CylDetector(detector.HoleRadius, detector.HoleDepth)
   	detout::CylDetector = CylDetector(detector.CryRadius, detector.CryLength)
   	Point(Height, pnt); #0.0 == SrcRadius && Point(pnt, 0.0)

   	if Height > 0.0							# the source as a whole out of the detector
      		return geoEff(detout, Point(Height, pnt), SrcRadius, SrcLength)

   	elseif Height + SrcLength < 0.0 		# the source as a whole inside of the detector
      		return 1.0 - geoEff(detin, Point(-(Height + SrcLength), pnt), SrcRadius, SrcLength)

   	elseif SrcLength > 0.0
      		res = (1.0 - geoEff(detin, Point(0.0), SrcRadius, -Height)) * -Height / SrcLength
      		res += geoEff(detout, Point(0.0), SrcRadius, Height + SrcLength) * (1.0 + Height / SrcLength)
      		return res

   	else
      		return 1.0 - geoEff(detin, Point(-Height, pnt), SrcRadius)

   	end #if
end #function

"""

	geoEff(detector::Detector, aPnt::Point, SrcRadius::Real = 0.0, SrcLength::Real = 0.0)::Float64

return the `geometrical efficiency` for a source (`point`, `disk` or `cylinder`) with
the detector `detector`.

## Arguments
*  `detector` can be any of the leaf detectors types (`CylDetector`, `BoreDetector`, `WellDetector`).
*  `aPNT`: a point represent the anchoring point of the source.
*  `SrcRadius`: Radius of the source.
*  `srcHeight`:  the height of an upright cylinder source.

## Throw
*  an `InValidGeometry` if the point location is invalide.
*  an `NotImplementedError` if source-to-detector geometry not supported yet.

!!! warning
    the point height of `aPnt` is measured differently for different detectors types.
    for the details, please refer to each detector entry.

!!! note
    *  if `SrcLength` equal to `zero`; the method return Geometrical Efficiency of a disc
       source of Radius = `SrcRadius` and center at the point `aPNT`.
    *  if both `SrcRadius` and `SrcLength` equal to `zero`;
       the method returns the Geometrical Efficiency of a point source at the anchoring point.

# Example

*   to obtain the efficiency of a `cylindrical` detector of crystal radius `2.0` cm for axial
    source cylinder of radius `1.0` cm and height `2.5` cm on the detector surface.

```jldoctest
julia> using GeoEfficiency

julia> geoEff(CylDetector(2.0), Point(0.0), 1.0, 2.5)
0.2923777934922748
```

*   to obtain the efficiency for a `bore-hole` detector of crystal radius of `2.0` and height of `3.0` with
    hole radius of `1.5` cm for axial source cylinder of radius `1.0` cm and height `2.5` cm starting from detector center.

```jldoctest
julia> using GeoEfficiency

julia> newDet = BoreDetector(2.0, 3.0, 1.5);

julia> geoEff(newDet, Point(0.0), 1.0, 2.5)
0.5678174038944723
```

*   to obtain the efficiency for a `well-type` detector of crystal radius of `2.0` cm and
    height `3.0` cm with hole radius of `1.5` cm and depth of `1.0` cm for axial source cylinder of
    radius `1.0` cm and height `2.5` cm at the hole surface.

```jldoctest
julia> using GeoEfficiency

julia> newDet = WellDetector(2.0, 3.0, 1.5, 1.0);

julia> geoEff(newDet, Point(0.0), 1.0, 2.5)
0.4669614527701105
```

"""
geoEff

"""

	geoEff(detector::Detector = Detector(), aSource::Tuple{Point, Real, Real} = source())::Float64

same as `geoEff(::Detector, ::Point, ::Real, ::Real)` but splatting the argument
in the Tuple `aSource`.

!!! note
    in the case of both no detector and no source Tuple or just the source Tuple is not supplied,
    it prompt the user to input a source (and detector) via the `console`.

"""
geoEff(detector::Detector = Detector(), aSource::Tuple{Point,Float64,Float64,} = source())::Float64 = geoEff(detector, aSource...)
geoEff(detector::Detector, aSource::Tuple{Point,Real,Real,})::Float64 = geoEff(detector, aSource...)


text:
"""
    Partitions a dataset.

    For each row in the dataset, check if it matches the question. If
    so, add it to 'true rows', otherwise, add it to 'false rows'.
"""
function partition(rows::Data, question::Question)::Tuple{Data, Data}
    true_rows, false_rows = Data(), Data()
    for row in rows
        if match(question, row)
            push!(true_rows, row)
        else
            push!(false_rows,row)
        end
    end
    true_rows, false_rows
end


text:
import LinearAlgebra: svd
import Base: getproperty

@explicit_intercepts svd Tuple{AbstractMatrix{<:Real}}

∇(::typeof(svd), ::Type{Arg{1}}, p, Y::SVD, Ȳ::NamedTuple{(:U,:S,:V)}, A::AbstractMatrix) =
    svd_rev(Y, Ȳ.U, Ȳ.S, Ȳ.V)

@explicit_intercepts getproperty Tuple{SVD, Symbol} [true, false]

function ∇(::typeof(getproperty), ::Type{Arg{1}}, p, y, ȳ, USV::SVD, x::Symbol)
    if x === :S
        return (U=zeroslike(USV.U), S=vec(ȳ), V=zeroslike(USV.V))
    elseif x === :U
        return (U=reshape(ȳ, size(USV.U)), S=zeroslike(USV.S), V=zeroslike(USV.V))
    elseif x === :V
        return (U=zeroslike(USV.U), S=zeroslike(USV.S), V=reshape(ȳ, size(USV.V)))
    elseif x === :Vt
        throw(ArgumentError("Vt is unsupported; use V and transpose the result"))
    else
        throw(ArgumentError("unrecognized property $x; expected U, S, or V"))
    end
end

function ∇(
    x̄::NamedTuple{(:U,:S,:V)},
    ::typeof(getproperty),
    ::Type{Arg{1}},
    p, y, ȳ,
    USV::SVD,
    x::Symbol,
)
    # This call does the validation that `x` is a recognized property
    x̄_update = ∇(getproperty, Arg{1}, p, y, ȳ, USV, x)
    if x === :S
        return (U=x̄.U, S=update!(x̄.S, x̄_update.S), V=x̄.V)
    elseif x === :U
        return (U=update!(x̄.U, x̄_update.U), S=x̄.S, V=x̄.V)
    elseif x === :V
        return (U=x̄.U, S=x̄.S, V=update!(x̄.V, x̄_update.V))
    end
end

# Iteration allows destructuring, e.g. U, S, V = svd(A)
# These definitions mirror those defined in the LinearAlgebra module, see
# https://github.com/JuliaLang/julia/blob/master/stdlib/LinearAlgebra/src/svd.jl#L20-L24
Base.iterate(usv::Branch{<:SVD}) = (usv.U, Val(:S))
Base.iterate(usv::Branch{<:SVD}, ::Val{:S}) = (usv.S, Val(:V))
Base.iterate(usv::Branch{<:SVD}, ::Val{:V}) = (usv.V, Val(:done))
Base.iterate(usv::Branch{<:SVD}, ::Val{:done}) = nothing

"""
    svd_rev(USV, Ū, S̄, V̄)

Compute the reverse mode sensitivities of the singular value decomposition (SVD). `USV` is
an `SVD` factorization object produced by a call to `svd`, and `Ū`, `S̄`, and `V̄` are the
respective sensitivities of the `U`, `S`, and `V` factors.
"""
function svd_rev(USV::SVD, Ū::AbstractMatrix, s̄::AbstractVector, V̄::AbstractMatrix)
    # Note: assuming a thin factorization, i.e. svd(A, full=false), which is the default
    U = USV.U
    s = USV.S
    V = USV.V
    Vt = USV.Vt

    k = length(s)
    T = eltype(s)
    F = T[i == j ? 1 : inv(@inbounds s[j]^2 - s[i]^2) for i = 1:k, j = 1:k]

    # We do a lot of matrix operations here, so we'll try to be memory-friendly and do
    # as many of the computations in-place as possible. Benchmarking shows that the in-
    # place functions here are significantly faster than their out-of-place, naively
    # implemented counterparts, and allocate no additional memory.
    Ut = U'
    FUᵀŪ = mulsubtrans!(Ut*Ū, F)  # F .* (UᵀŪ - ŪᵀU)
    FVᵀV̄ = mulsubtrans!(Vt*V̄, F)  # F .* (VᵀV̄ - V̄ᵀV)
    ImUUᵀ = eyesubx!(U*Ut)        # I - UUᵀ
    ImVVᵀ = eyesubx!(V*Vt)        # I - VVᵀ

    S = Diagonal(s)
    S̄ = Diagonal(s̄)

    Ā = add!(U*FUᵀŪ*S, ImUUᵀ*(Ū/S))*Vt
    add!(Ā, U*S̄*Vt)
    add!(Ā, U*add!(S*FVᵀV̄*Vt, (S\V̄')*ImVVᵀ))

    return Ā
end

"""
    mulsubtrans!(X::AbstractMatrix, F::AbstractMatrix)

Compute `F .* (X - X')`, overwriting `X` in the process.

!!! note
    This is an internal function that does no argument checking; the matrices passed to
    this function are square with matching dimensions by construction.
"""
function mulsubtrans!(X::AbstractMatrix{T}, F::AbstractMatrix{T}) where T<:Real
    k = size(X, 1)
    @inbounds for j = 1:k, i = 1:j  # Iterate the upper triangle
        if i == j
            X[i,i] = zero(T)
        else
            X[i,j], X[j,i] = F[i,j] * (X[i,j] - X[j,i]), F[j,i] * (X[j,i] - X[i,j])
        end
    end
    X
end

"""
    eyesubx!(X::AbstractMatrix)

Compute `I - X`, overwriting `X` in the process.
"""
function eyesubx!(X::AbstractMatrix{T}) where T<:Real
    n, m = size(X)
    @inbounds for j = 1:m, i = 1:n
        X[i,j] = (i == j) - X[i,j]
    end
    X
end

"""
    add!(X::AbstractMatrix, Y::AbstractMatrix)

Compute `X + Y`, overwriting X in the process.

!!! note
    This is an internal function that does no argument checking; the matrices passed to
    this function are square with matching dimensions by construction.
"""
function add!(X::AbstractMatrix{T}, Y::AbstractMatrix{T}) where T<:Real
    @inbounds for i = eachindex(X, Y)
        X[i] += Y[i]
    end
    X
end


text:
# Julia wrapper for header: cpl_virtualmem.h
# Automatically generated using Clang.jl

"""
    CPLGetPageSize(void) -> size_t

Return the size of a page of virtual memory.

### Returns
the page size.
"""
function cplgetpagesize()
    aftercare(ccall((:CPLGetPageSize, libgdal), Csize_t, ()))
end

"""
    CPLVirtualMemNew(size_t nSize,
                     size_t nCacheSize,
                     size_t nPageSizeHint,
                     int bSingleThreadUsage,
                     CPLVirtualMemAccessMode eAccessMode,
                     CPLVirtualMemCachePageCbk pfnCachePage,
                     CPLVirtualMemUnCachePageCbk pfnUnCachePage,
                     CPLVirtualMemFreeUserData pfnFreeUserData,
                     void * pCbkUserData) -> CPLVirtualMem *

Create a new virtual memory mapping.

### Parameters
* **nSize**: size in bytes of the virtual memory mapping.
* **nCacheSize**: size in bytes of the maximum memory that will be really allocated (must ideally fit into RAM).
* **nPageSizeHint**: hint for the page size. Must be a multiple of the system page size, returned by CPLGetPageSize(). Minimum value is generally 4096. Might be set to 0 to let the function determine a default page size.
* **bSingleThreadUsage**: set to TRUE if there will be no concurrent threads that will access the virtual memory mapping. This can optimize performance a bit.
* **eAccessMode**: permission to use for the virtual memory mapping.
* **pfnCachePage**: callback triggered when a still unmapped page of virtual memory is accessed. The callback has the responsibility of filling the page with relevant values.
* **pfnUnCachePage**: callback triggered when a dirty mapped page is going to be freed (saturation of cache, or termination of the virtual memory mapping). Might be NULL.
* **pfnFreeUserData**: callback that can be used to free pCbkUserData. Might be NULL
* **pCbkUserData**: user data passed to pfnCachePage and pfnUnCachePage.

### Returns
a virtual memory object that must be freed by CPLVirtualMemFree(), or NULL in case of failure.
"""
function cplvirtualmemnew(nSize, nCacheSize, nPageSizeHint, bSingleThreadUsage, eAccessMode, pfnCachePage, pfnUnCachePage, pfnFreeUserData, pCbkUserData)
    aftercare(ccall((:CPLVirtualMemNew, libgdal), Ptr{CPLVirtualMem}, (Csize_t, Csize_t, Csize_t, Cint, CPLVirtualMemAccessMode, CPLVirtualMemCachePageCbk, CPLVirtualMemUnCachePageCbk, CPLVirtualMemFreeUserData, Ptr{Cvoid}), nSize, nCacheSize, nPageSizeHint, bSingleThreadUsage, eAccessMode, pfnCachePage, pfnUnCachePage, pfnFreeUserData, pCbkUserData))
end

"""
    CPLIsVirtualMemFileMapAvailable(void) -> int

Return if virtual memory mapping of a file is available.

### Returns
TRUE if virtual memory mapping of a file is available.
"""
function cplisvirtualmemfilemapavailable()
    aftercare(ccall((:CPLIsVirtualMemFileMapAvailable, libgdal), Cint, ()))
end

"""
    CPLVirtualMemFileMapNew(VSILFILE * fp,
                            vsi_l_offset nOffset,
                            vsi_l_offset nLength,
                            CPLVirtualMemAccessMode eAccessMode,
                            CPLVirtualMemFreeUserData pfnFreeUserData,
                            void * pCbkUserData) -> CPLVirtualMem *

Create a new virtual memory mapping from a file.

### Parameters
* **fp**: Virtual file handle.
* **nOffset**: Offset in the file to start the mapping from.
* **nLength**: Length of the portion of the file to map into memory.
* **eAccessMode**: Permission to use for the virtual memory mapping. This must be consistent with how the file has been opened.
* **pfnFreeUserData**: callback that is called when the object is destroyed.
* **pCbkUserData**: user data passed to pfnFreeUserData.

### Returns
a virtual memory object that must be freed by CPLVirtualMemFree(), or NULL in case of failure.
"""
function cplvirtualmemfilemapnew(fp, nOffset, nLength, eAccessMode, pfnFreeUserData, pCbkUserData)
    aftercare(ccall((:CPLVirtualMemFileMapNew, libgdal), Ptr{CPLVirtualMem}, (Ptr{VSILFILE}, vsi_l_offset, vsi_l_offset, CPLVirtualMemAccessMode, CPLVirtualMemFreeUserData, Ptr{Cvoid}), fp, nOffset, nLength, eAccessMode, pfnFreeUserData, pCbkUserData))
end

"""
    CPLVirtualMemDerivedNew(CPLVirtualMem * pVMemBase,
                            vsi_l_offset nOffset,
                            vsi_l_offset nSize,
                            CPLVirtualMemFreeUserData pfnFreeUserData,
                            void * pCbkUserData) -> CPLVirtualMem *

Create a new virtual memory mapping derived from an other virtual memory mapping.

### Parameters
* **pVMemBase**: Base virtual memory mapping
* **nOffset**: Offset in the base virtual memory mapping from which to start the new mapping.
* **nSize**: Size of the base virtual memory mapping to expose in the the new mapping.
* **pfnFreeUserData**: callback that is called when the object is destroyed.
* **pCbkUserData**: user data passed to pfnFreeUserData.

### Returns
a virtual memory object that must be freed by CPLVirtualMemFree(), or NULL in case of failure.
"""
function cplvirtualmemderivednew(pVMemBase, nOffset, nSize, pfnFreeUserData, pCbkUserData)
    aftercare(ccall((:CPLVirtualMemDerivedNew, libgdal), Ptr{CPLVirtualMem}, (Ptr{CPLVirtualMem}, vsi_l_offset, vsi_l_offset, CPLVirtualMemFreeUserData, Ptr{Cvoid}), pVMemBase, nOffset, nSize, pfnFreeUserData, pCbkUserData))
end

"""
    CPLVirtualMemFree(CPLVirtualMem * ctxt) -> void

Free a virtual memory mapping.

### Parameters
* **ctxt**: context returned by CPLVirtualMemNew().
"""
function cplvirtualmemfree(ctxt)
    aftercare(ccall((:CPLVirtualMemFree, libgdal), Cvoid, (Ptr{CPLVirtualMem},), ctxt))
end

"""
    CPLVirtualMemGetAddr(CPLVirtualMem * ctxt) -> void *

Return the pointer to the start of a virtual memory mapping.

### Parameters
* **ctxt**: context returned by CPLVirtualMemNew().

### Returns
the pointer to the start of a virtual memory mapping.
"""
function cplvirtualmemgetaddr(ctxt)
    aftercare(ccall((:CPLVirtualMemGetAddr, libgdal), Ptr{Cvoid}, (Ptr{CPLVirtualMem},), ctxt))
end

"""
    CPLVirtualMemGetSize(CPLVirtualMem * ctxt) -> size_t

Return the size of the virtual memory mapping.

### Parameters
* **ctxt**: context returned by CPLVirtualMemNew().

### Returns
the size of the virtual memory mapping.
"""
function cplvirtualmemgetsize(ctxt)
    aftercare(ccall((:CPLVirtualMemGetSize, libgdal), Csize_t, (Ptr{CPLVirtualMem},), ctxt))
end

"""
    CPLVirtualMemIsFileMapping(CPLVirtualMem * ctxt) -> int

Return if the virtual memory mapping is a direct file mapping.

### Parameters
* **ctxt**: context returned by CPLVirtualMemNew().

### Returns
TRUE if the virtual memory mapping is a direct file mapping.
"""
function cplvirtualmemisfilemapping(ctxt)
    aftercare(ccall((:CPLVirtualMemIsFileMapping, libgdal), Cint, (Ptr{CPLVirtualMem},), ctxt))
end

"""
    CPLVirtualMemGetAccessMode(CPLVirtualMem * ctxt) -> CPLVirtualMemAccessMode

Return the access mode of the virtual memory mapping.

### Parameters
* **ctxt**: context returned by CPLVirtualMemNew().

### Returns
the access mode of the virtual memory mapping.
"""
function cplvirtualmemgetaccessmode(ctxt)
    aftercare(ccall((:CPLVirtualMemGetAccessMode, libgdal), CPLVirtualMemAccessMode, (Ptr{CPLVirtualMem},), ctxt))
end

"""
    CPLVirtualMemGetPageSize(CPLVirtualMem * ctxt) -> size_t

Return the page size associated to a virtual memory mapping.

### Parameters
* **ctxt**: context returned by CPLVirtualMemNew().

### Returns
the page size
"""
function cplvirtualmemgetpagesize(ctxt)
    aftercare(ccall((:CPLVirtualMemGetPageSize, libgdal), Csize_t, (Ptr{CPLVirtualMem},), ctxt))
end

"""
    CPLVirtualMemIsAccessThreadSafe(CPLVirtualMem * ctxt) -> int

Return TRUE if this memory mapping can be accessed safely from concurrent threads.

### Parameters
* **ctxt**: context returned by CPLVirtualMemNew().

### Returns
TRUE if this memory mapping can be accessed safely from concurrent threads.
"""
function cplvirtualmemisaccessthreadsafe(ctxt)
    aftercare(ccall((:CPLVirtualMemIsAccessThreadSafe, libgdal), Cint, (Ptr{CPLVirtualMem},), ctxt))
end

"""
    CPLVirtualMemDeclareThread(CPLVirtualMem * ctxt) -> void

Declare that a thread will access a virtual memory mapping.

### Parameters
* **ctxt**: context returned by CPLVirtualMemNew().
"""
function cplvirtualmemdeclarethread(ctxt)
    aftercare(ccall((:CPLVirtualMemDeclareThread, libgdal), Cvoid, (Ptr{CPLVirtualMem},), ctxt))
end

"""
    CPLVirtualMemUnDeclareThread(CPLVirtualMem * ctxt) -> void

Declare that a thread will stop accessing a virtual memory mapping.

### Parameters
* **ctxt**: context returned by CPLVirtualMemNew().
"""
function cplvirtualmemundeclarethread(ctxt)
    aftercare(ccall((:CPLVirtualMemUnDeclareThread, libgdal), Cvoid, (Ptr{CPLVirtualMem},), ctxt))
end

"""
    CPLVirtualMemPin(CPLVirtualMem * ctxt,
                     void * pAddr,
                     size_t nSize,
                     int bWriteOp) -> void

Make sure that a region of virtual memory will be realized.

### Parameters
* **ctxt**: context returned by CPLVirtualMemNew().
* **pAddr**: the memory region to pin.
* **nSize**: the size of the memory region.
* **bWriteOp**: set to TRUE if the memory are will be accessed in write mode.
"""
function cplvirtualmempin(ctxt, pAddr, nSize, bWriteOp)
    aftercare(ccall((:CPLVirtualMemPin, libgdal), Cvoid, (Ptr{CPLVirtualMem}, Ptr{Cvoid}, Csize_t, Cint), ctxt, pAddr, nSize, bWriteOp))
end

"""
    CPLVirtualMemManagerTerminate(void) -> void

Cleanup any resource and handlers related to virtual memory.
"""
function cplvirtualmemmanagerterminate()
    aftercare(ccall((:CPLVirtualMemManagerTerminate, libgdal), Cvoid, ()))
end


text:
module RWUtilities

using FFTW
import Random
using LinearAlgebra

import Printf

#import HCubature

import GSL
import SpecialFunctions

import Interpolations
#import StaticArrays
import StatsFuns

import Test

include("declarations.jl")

include("matrix.jl")
include("approximations.jl")
include("geometric_transforms.jl")
include("signals.jl")
include("boundary.jl")
include("patterns.jl")
include("bounds.jl")
include("combinatorics.jl")
include("template_signals.jl")
include("random_matrices.jl")
include("data_structures.jl")
include("simplex.jl")
include("function_interpolation.jl")
include("numerical.jl")

export  isnumericallyclose, # bounds.jl
        isnumericallyin,
        makenotclose!,

        # combinatorics
        coprime,
        EuclideanAlgotihm,

        # patterns.jl
        extractselectcols,
        finddiscrepancy,
        forcesymmetric,
        forcenegativesymmetric,
        coordinatetolattice1D,
        lattice1Dtocoordinate,
        ranges2collection,

        # boundary.jl
        extractframecenters,
        padzeros,
        evalzeropaddedsignal,

        # signals.jl
        rad2freq,
        freq2rad,
        computeDTFTviaformula,
        getdiracdeltaseq,
        spectralinversion,
        upsample,
        downsample,
        linearconv,
        forcesamesizeafterlinearconv,
        fastlinearconv,
        circconv,

        # random_matrices.jl
        generaterandomsym,
        generaterandomskew,
        generaterandomortho,
        generateStiefel,
        generaterandomhankel,
        generateindefmatrix,
        generateposdefmatrix,

        # template_signals.jl
        rectfunc,
        raisedcosinefunction,
        cosinetransitionlowpassfunction,
        cosinetransitionhighpassfunction,
        cosinetransitionbandpassfunction,

        # geometric_transforms.jl
        shiftthenstretch,
        shiftthencompress,
        evalstretchedsignal,
        convertcompactdomain,

        # approximations.jl
        smoothmin,
        smoothmax,
        smoothextrema,

        # matrix.jl
        forcesymmetric,
        naivesqrtpsdmatrix,
        array2matrix,
        makeblockdiagonalmatrix,
        invertblockmatrix,
        blockmatrixmatrixproduct,
        blockmatrixvectorproduct,

        # data_structures.jl
        packbcW,
        unpackbcW,
        parseCholeskyType,
        packCholeskyType,

        # simplex.jl
        simplexinteriortoEuclidean,
        Euclideantosimplex,

        # function_interpolation
        setupcubicitp,

        # numerical
        logsumexp
end

max_stars_repo_name: RoyCCWang/RWUtilities

max_issues_repo_name: RoyCCWang/RWUtilities

max_forks_repo_name: RoyCCWang/RWUtilities


text:
using Distributions
using MonteCarloMeasurements

export importanceSample
@inline function importanceSample(p::JointDistribution, q::JointDistribution, _data)
    return _importanceSample(p.model, p.args, q.model, q.args, _data)
end

@gg function _importanceSample(p::Model, _pargs, q::Model, _qargs, _data)
    p = type2model(p)
    q = type2model(q)

    sourceImportanceSample()(p,q) |> loadvals(_qargs, _data) |> loadvals(_pargs, NamedTuple())
end

export sourceImportanceSample
function sourceImportanceSample()
    function(p::Model,q::Model)
        p = canonical(p)
        q = canonical(q)
        m = merge(p,q)

        function proc(m, st::Sample)
            if hasproperty(p.dists, st.x)
                pdist = getproperty(p.dists, st.x)
                qdist = st.rhs
                @gensym ℓx
                result = @q begin
                    $ℓx = importanceSample($pdist, $qdist, _data)
                    _ℓ += $ℓx.ℓ
                    $(st.x) = $ℓx.val
                end
                return flatten(result)
            else return :($(st.x) = $(st.rhs))
            end
            return flatten(result)
        end
        proc(m, st::Assign)     = :($(st.x) = $(st.rhs))
        proc(m, st::Return)  = :(return $(st.rhs))
        proc(m, st::LineNumber) = nothing

        body = buildSource(m, proc) |> flatten

        kwargs = freeVariables(q) ∪ arguments(p)
        kwargsExpr = Expr(:tuple,kwargs...)

        stochExpr = begin
            vals = map(stochastic(m)) do x Expr(:(=), x,x) end
            Expr(:tuple, vals...)
        end

        wrap(kernel) = @q begin
            _ℓ = 0.0
            $body
            return Weighted(_ℓ, $stochExpr)
        end

        buildSource(m, proc, wrap) |> flatten
    end
end

@inline function importanceSample(p, q, _data)
    x = merge(rand(q), _data)
    ℓ = logpdf(p,x) - logpdf(q,x)
    Weighted(ℓ,x)
end

# export sourceParticleImportance
# function sourceParticleImportance(p,q;ℓ=:ℓ)
#     p = canonical(p)
#     q = canonical(q)
#     @gensym N
#     # This determines how to initialize a Particle for a given expression
#     vars(expr) = (bound(p) ∪ bound(q) ∪ stochastic(p) ∪ stochastic(q)) ∩ variables(expr)

#     procp(p, st::Follows) = :($ℓ += logpdf($(st.rhs), $(st.x)))
#     procp(p, st::Let)     = convert(Expr, st)
#     procp(p, st::Return)  = nothing
#     procp(p, st::LineNumber) = convert(Expr, st)

#     function procq(q, st::Follows)
#         if isempty(vars(st.rhs))
#             @q begin
#                 $(st.x) = Particles($N, $(st.rhs))
#                 $ℓ -= logpdf($(st.rhs), $(st.x))
#             end
#         else
#             @q begin
#                 $(st.x) = rand($(st.rhs))
#                 $ℓ -= logpdf($(st.rhs), $(st.x))
#             end
#         end
#     end
#     procq(q, st::Let)     = convert(Expr, st)
#     procq(q, st::Return)  = nothing
#     procq(q, st::LineNumber) = convert(Expr, st)

#     pbody = buildSource(p, procp) |> striplines
#     qbody = buildSource(q, procq) |> striplines

#     kwargs = freeVariables(q) ∪ arguments(p)
#     kwargsExpr = Expr(:tuple,kwargs...)

#     stochExpr = begin
#         vals = map(stochastic(q)) do x Expr(:(=), x,x) end
#         Expr(:tuple, vals...)
#     end

#     @gensym particleImportance
#     result = @q function $particleImportance($N, pars)
#         @unpack $kwargsExpr = pars
#         $ℓ = 0.0 * Particles($N, Uniform())
#         $qbody
#         $pbody
#         return ($ℓ, $stochExpr)
#     end

#     flatten(result)
# end

#     @gensym rand

#     flatten(@q (
#         function $rand(args...;kwargs...)
#             @unpack $argsExpr = kwargs
#             # kwargs = Dict(kwargs)
#             $body
#             $stochExpr
#         end
#     ))

# end


text:
struct Path
    density::Vector{Float64}
    baseline::Vector{Float64}
end

Path(density::Number, baseline::Number) = Path([density],[baseline])

Base.iterate(p::Path, state=1) = state > length(p.density) ? nothing : ( (p.density[state], p.baseline[state]),  state+1)

Base.length(p::Path) = length(p.density)

"""
$(SIGNATURES)

Create modified oscillation parameters for neutrino propagation through matter

# Arguments
- `P`: Vacuum PMNS Matrix
- `H`: Vacuum Hamiltonian
- `density`: Matter density [g*cm^-3]
- `energy`: Neutrino energy [GeV]
- `zoa`: Proton nucleon ratio (Z/A)
- `anti`: Is anti neutrino
"""
function MatterOscillationMatrices(U, H, energy, density; zoa=0.5, anti=false)
    H_eff = convert(Array{ComplexF64}, U * Diagonal{Complex}(H) * adjoint(U))
    MatterOscillationMatrices(H_eff, energy, density; zoa=zoa, anti=anti)
end

"""
$(SIGNATURES)

Create modified oscillation parameters for neutrino propagation through matter

# Arguments
- `H_eff`: Effective Matter Hamiltonian
- `density`: Matter density [g*cm^-3]
- `energy`: Neutrino energy [GeV]
- `zoa`: Proton nucleon ratio (Z/A)
- `anti`: Is anti neutrino
"""
function MatterOscillationMatrices(H_eff, energy, density; zoa=0.5, anti=false)
    A = sqrt(2) * G_F * N_A * density
    if anti
        H_eff[1,1] -= A * zoa * 2 * energy * 1e9
    else
        H_eff[1,1] += A * zoa * 2 * energy * 1e9
    end
    # Subtract Vz (= - A*Nn*E*1e9) for any sterile flavours
    if size(H_eff)[1] > 3
        for i in 4:size(H_eff)[1]
            if anti
                H_eff[i,i] -= A * (1 - zoa) * energy * 1e9
            else
                H_eff[i,i] += A * (1 - zoa) * energy * 1e9
            end
        end
    end
    tmp = eigen(H_eff)
    return tmp.vectors, tmp.values
end

"""
$(SIGNATURES)

Create modified oscillation parameters for neutrino propagation through matter

# Arguments
- `osc_vacuum::OscillationParameters`: Oscillation parameters in vacuum
- `energy`: Neutrino energy [GeV]
- `density`: Matter density in g*cm^-3
- `zoa`: Proton nucleon ratio (Z/A)
- `anti`: Is anti neutrino

"""
function MatterOscillationMatrices(osc_vacuum::OscillationParameters, energy, density; zoa=0.5, anti=false)
    H_vacuum = Diagonal(Hamiltonian(osc_vacuum))
    U_vacuum = PMNSMatrix(osc_vacuum; anti=anti)
    H_eff = convert(Array{ComplexF64}, U_vacuum * Diagonal{ComplexF64}(H_vacuum) * adjoint(U_vacuum))
    return MatterOscillationMatrices(H_eff, energy, density; zoa=zoa, anti=anti)
end

"""
$(SIGNATURES)

# Arguments
- `U`: Vacuum PMNS Matrix
- `H`: Vacuum Hamiltonian
- `energy`: Neutrino energy [GeV]
- `path::Vector{Path}`: Neutrino path
- `zoa`: Proton nucleon ratio (Z/A)
- `anti`: Is anti neutrino
"""
function oscprob(U, H, energy::Vector{T}, path::Vector{Path}; zoa=0.5, anti=false) where {T <: Real}
    energy = convert.(Float64, energy)
    if anti
        H_eff = conj.(U) * Diagonal{ComplexF64}(H) * adjoint(conj.(U))
    else
        H_eff = U * Diagonal{ComplexF64}(H) * adjoint(U)
    end
    A = zeros(ComplexF64, length(energy), length(path), size(U)...)
    cache_size = length(energy) * sum(map(x->length(x.density), path))
    lru = LRU{Tuple{Float64, Float64},
              Tuple{Array{ComplexF64,2}, Vector{ComplexF64}}}(maxsize=cache_size)
    for k in 1:length(energy)
        @inbounds E = energy[k]
        for (l, p) in enumerate(path)
            tmp = Matrix{ComplexF64}(1I, size(U))
            for (m,b) in enumerate(p.baseline)
                @inbounds ρ = p.density[m]
                U_mat, H_mat = get!(lru, (E, ρ)) do
                    MatterOscillationMatrices(copy(H_eff), E, ρ; zoa=zoa, anti=anti)
                end
                tmp *= Neurthino._oscprobampl(U_mat, H_mat, E, b)
            end
            @inbounds A[k, l,  :, :] = tmp
        end
    end
    P = map(x -> abs.(x) .^ 2, A)
    flavrange = _make_flavour_range(first(size(U)))
    AxisArray(P; Energy=energy, Path=path, InitFlav=flavrange, FinalFlav=flavrange)
end

const oscprob(U, H, energy::T, path::Vector{Path}; zoa=0.5, anti=false) where {T <: Real} = oscprob(U, H, [energy], path; zoa=zoa, anti=anti)

const oscprob(U, H, energy, path::Path; zoa=0.5, anti=false) = oscprob(U, H, energy, [path]; zoa=zoa, anti=anti)

"""
$(SIGNATURES)

# Arguments
- `osc_vacuum::OscillationParameters`: Vacuum oscillation parameters
- `energy`: Neutrino energy [GeV]
- `path`: Neutrino path
- `zoa`: Proton nucleon ratio (Z/A)
- `anti`: Is anti neutrino
"""
function oscprob(osc_vacuum::OscillationParameters, energy, path::Union{Path, Vector{Path}}; zoa=0.5, anti=false)
    # TODO: attach U_vac and H_vac to the oscillation parameters, so that it's
    # only calculated once and invalidated when any of the oscillation parameters
    # are changed
    U_vac = PMNSMatrix(osc_vacuum; anti=anti)
    H_vac = Hamiltonian(osc_vacuum)
    oscprob(U_vac, H_vac, energy, path; zoa=zoa, anti=anti)
end


text:
# Autogenerated wrapper script for FFMPEG_jll for x86_64-w64-mingw32
export ffmpeg, ffprobe, libavcodec, libavdevice, libavfilter, libavformat, libavresample, libavutil, libpostproc, libswresample, libswscale

using libass_jll
using libfdk_aac_jll
using FriBidi_jll
using FreeType2_jll
using LAME_jll
using libvorbis_jll
using Ogg_jll
using x264_jll
using x265_jll
using Bzip2_jll
using Zlib_jll
using OpenSSL_jll
using Opus_jll
JLLWrappers.@generate_wrapper_header("FFMPEG")
JLLWrappers.@declare_executable_product(ffmpeg)
JLLWrappers.@declare_executable_product(ffprobe)
JLLWrappers.@declare_library_product(libavcodec, "avcodec-58.dll")
JLLWrappers.@declare_library_product(libavdevice, "avdevice-58.dll")
JLLWrappers.@declare_library_product(libavfilter, "avfilter-7.dll")
JLLWrappers.@declare_library_product(libavformat, "avformat-58.dll")
JLLWrappers.@declare_library_product(libavresample, "avresample-4.dll")
JLLWrappers.@declare_library_product(libavutil, "avutil-56.dll")
JLLWrappers.@declare_library_product(libpostproc, "postproc-55.dll")
JLLWrappers.@declare_library_product(libswresample, "swresample-3.dll")
JLLWrappers.@declare_library_product(libswscale, "swscale-5.dll")
function __init__()
    JLLWrappers.@generate_init_header(libass_jll, libfdk_aac_jll, FriBidi_jll, FreeType2_jll, LAME_jll, libvorbis_jll, Ogg_jll, x264_jll, x265_jll, Bzip2_jll, Zlib_jll, OpenSSL_jll, Opus_jll)
    JLLWrappers.@init_executable_product(
        ffmpeg,
        "bin\\ffmpeg.exe",
    )

    JLLWrappers.@init_executable_product(
        ffprobe,
        "bin\\ffprobe.exe",
    )

    JLLWrappers.@init_library_product(
        libavcodec,
        "bin\\avcodec-58.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libavdevice,
        "bin\\avdevice-58.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libavfilter,
        "bin\\avfilter-7.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libavformat,
        "bin\\avformat-58.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libavresample,
        "bin\\avresample-4.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libavutil,
        "bin\\avutil-56.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libpostproc,
        "bin\\postproc-55.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libswresample,
        "bin\\swresample-3.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libswscale,
        "bin\\swscale-5.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
using GraphKernels
using Test

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


text:
type Timematr{T} <: AbstractTimematr
    vals::DataFrame
    idx::Array{T, 1}

    function Timematr(vals::DataFrame, idx::Array{T, 1})
        chkIdx(idx)
        chkNum(vals)
        if(size(vals, 1) != length(idx))
            if (length(idx) == 0) | (size(vals, 1) == 0)
                df = convert(DataFrame, rand(2, 2))
                return new(df[1:2, [false, false]], Array{T, 1}[])
            end
            error(length(idx), " idx entries, but ", size(vals, 1), " rows of data")
        end
        return new(vals, idx)
    end
end

function Timematr{T}(vals::DataFrame, idx::Array{T, 1})
    return Timematr{T}(vals, idx)
end

#############
## getVars ##
#############

function getVars(tm::Timematr, mapF::Function, crit::Function)
    ## map each column to a value, and apply some condition on value
    ## e.g.: find variables with minimum value less than -30
    ncol = size(tm, 2)
    nrow = size(tm, 1)
    vals = TimeData.core(tm)

    values = ones(ncol)
    for ii=1:ncol
        values[ii] = mapF(vals[:, ii])
    end

    logics = crit(values)

    df = DataFrame(names=names(tm)[logics], values=values[logics])
    ## names!(df, names(tm)[logics])
    return df
end


text:
# module HitTest

# using Gtk, GtkUtilities, ..Graphics, Colors
# import Compose
# import ..Immerse: nearest, hitcenter, Figure, render_backend, find_object

# export
#     hit,
#     circle_center

# import .DisplayGadfly: _hit_data

"""
`hit((fig,tag), cb)` turns on hit-testing for the plot element tagged with `tag` in figure `fig`. When the user clicks on the element, the callback function `cb` will be executed.

The callback should have the form

```jl
    function my_callback(figtag, index, event, distance)
        if distance < 2
            # We clicked close enough to the object to "count"
            # implement the action here
        end
    end
```

`figtag = (figno,tag)` identifies the figure and object, and `index`
represents the specific item clicked on.  `index` might be a single `Int`
(e.g., the 37th point in the plot), or for compound objects might be a
more complex object.  `event` contains Gtk's information about the
click event; this records the `x,y` position as well as any modifier
keys.  `distance` measures how much the user "missed" the object, and
is measured in screen-pixels.

Here are a couple of examples:
```jl
    fig = display(plot(x=1:10,y=rand(10),Geom.point(tag=:dots)))
    hit((fig,:dots), (figtag, index, xy, dist) -> if dist < 1 println("You clicked on dot ", index) end)
```
In this case, `index` will be an `Int`.

```jl
    x = rand(6)
    y = rand(6)
    label = [1,1,1,1,2,2]
    df = DataFrame(Any[x,y,label], [:x,:y,:label])
    fig = display(plot(df, x=:x, y=:y, color=:label, Geom.line(tag=:lines)))

    hit((fig,:lines), (figtag, index, xy, dist) -> begin
        if dist < 2
            println("You clicked on line ", index[1], " in segment ", index[2], " at \$(round(Int,100*index[3]))% along the segment")
        end
    end)
```

By default, `hit` uses mouse button 1; you can customize the action with `hit(figtag, cb, action)` where the default is `action = :button1press`.  The choices correspond to the fields in Gtk's `MouseHandler`.

You can suspend hit-testing for an object with `hit((fig,:tag), false)`
and re-enable it with `hit((fig,:tag), true)`.
"""
function hit(figtag::Tuple{Int,Symbol}, cb, action = :button1press)
    figno, tag = figtag
    fig = Figure(figno)
    if !haskey(_hit_data, fig)
        _hit_data[fig] = Dict{Symbol,Any}()
        c = fig.canvas
        push!((c.mouse, action),
              Gtk.@guarded (widget, event) -> hitcb(figno, fig, event))
    end
    _hit_data[fig][tag] = (true, cb)   # starts in "on" position
    nothing
end

function hit(figtag::Tuple{Int,Symbol}, state::Bool)
    figno, tag = figtag
    dct = _hit_data[Figure(figno)]
    olddata = dct[tag]
    dct[tag] = (state, olddata[2])
    nothing
end

# callback function for hit-testing
function hitcb(figno, f, event)
    c = f.canvas
    x, y = event.x, event.y
    hitables = get(_hit_data, f, nothing)
    hitables == nothing && return
    coords = GtkUtilities.guidata[c, :coords]
    mindist = Inf
    minindex = 0   # not type-stable, unfortunately
    obj = Compose.empty_tag
    backend = render_backend(c)
    local objcb
    for (tag, data) in hitables
        state, cb = data
        !state && continue
        # Find the object in the rendered figure
        form = find_object(f.cc, tag)
        dist, index = nearest(backend, coords[tag], form, x, y)
        if dist < mindist
            mindist = dist
            obj = tag
            objcb = cb
            minindex = index
        end
    end
    if obj != Compose.empty_tag
        objcb((figno,obj), minindex, event, mindist)
    end
    nothing
end

# A good callback function for testing
#    hit(figtag, (figtag, index, event, dist) -> circle_center(figtag, index))
function circle_center(figtag::Tuple{Int,Symbol}, index; color=RGB{U8}(1,0,0))
    figno, tag = figtag
    f = Figure(figno)
    c = f.canvas
    coords = GtkUtilities.guidata[c, :coords][tag]
    form = find_object(f.cc, tag)
    backend = render_backend(c)
    x, y = hitcenter(backend, coords, form, index)
    ctx = getgc(c)
    set_source(ctx, color)
    set_line_width(ctx, 2)
    arc(ctx, x, y, 5, 0, 2pi)
    stroke(ctx)
    Gtk.reveal(c)
    nothing
end

# end  # module


text:
#ДАНО: На ограниченном внешней прямоугольной рамкой поле имеется ровно одна внутренняя перегородка в форме прямоугольника.
#Робот - в произвольной клетке поля между внешней и внутренней перегородками.
#РЕЗУЛЬТАТ: Робот - в исходном положении и по всему периметру внутренней перегородки поставлены маркеры.
function move_to_startplace!(r::Robot)
    x = 0
    y = 0
    while !(isborder(r, HorizonSide(2)) && isborder(r, HorizonSide(3)))
        while !isborder(r, HorizonSide(2))
            move!(r, HorizonSide(2))
            y+=1
        end
        while !isborder(r, HorizonSide(3))
            move!(r, HorizonSide(3))
            x+=1
        end
    end
    return (x, y)
end

function move_to_beginplace!(r::Robot, X, Y)
    x = X
    y = Y

    while x != 0 || y != 0
        while x != 0
            if x > 0
                if isborder(r, HorizonSide(1))
                    move!(r, HorizonSide(0))
                    y-=1
                end
                if !isborder(r, HorizonSide(1))
                    move!(r, HorizonSide(1))
                    x -= 1
                end
            elseif x < 0
                if isborder(r, HorizonSide(3))
                    move!(r, HorizonSide(0))
                    y-=1
                end
                if !isborder(r, HorizonSide(3))
                    move!(r, HorizonSide(3))
                    x += 1
                end
            end
        end
        while y != 0

            if y > 0
                if isborder(r, HorizonSide(0))
                    move!(r, HorizonSide(1))
                    x-=1
                end
                if !isborder(r, HorizonSide(0))
                    move!(r, HorizonSide(0))
                    y -= 1
                end
            elseif y < 0
                if isborder(r, HorizonSide(2))
                    move!(r, HorizonSide(1))
                    x-=1
                end
                if !isborder(r, HorizonSide(2))
                    move!(r, HorizonSide(2))
                    y += 1
                end
            end
        end
    end
end

function search_object!(r::Robot)
    x = 0
    y = 0
    i = 2
    check = false
    while !isborder(r, HorizonSide(0))
        move!(r, HorizonSide(0))
        y+=1
    end
    while !isborder(r, HorizonSide(1))
        move!(r, HorizonSide(1))
        x+=1
    end
    while !check
        if i % 2 == 0
            for _ in 1:(y-1)
                move!(r, HorizonSide(i))
                if isborder(r, HorizonSide((i+1)%4)) || isborder(r, HorizonSide(i))
                    check = true
                    return i
                end
            end
            y-=1
        end
        if i % 2 != 0
            for _ in 1:(x-1)
                move!(r, HorizonSide(i))
                if isborder(r, HorizonSide((i+1)%4)) || isborder(r, HorizonSide(i))
                    check = true
                    return i
                end
            end
            x-=1
        end
        i = (i + 1)%4
    end
    return ErrorException
end

function perimetr_around_object!(r::Robot, side)
    check = false

    while !check
        if ismarker(r)
            break
        end
        while isborder(r, HorizonSide((side+1)%4))
            putmarker!(r)
            move!(r, HorizonSide(side))
        end
        putmarker!(r)
        side = (side+1)%4
        move!(r, HorizonSide(side))
    end
end

function okolo_obyekta(r::Robot)
    x, y = move_to_startplace!(r)
    side = search_object!(r)
    perimetr_around_object!(r, side)
    move_to_startplace!(r)
    move_to_beginplace!(r, x, y)
end

max_stars_repo_name: listova-nastya/mirea-progs

max_issues_repo_name: listova-nastya/mirea-progs

max_forks_repo_name: listova-nastya/mirea-progs


text:
# include(joinpath(dirname(dirname(@__FILE__)), "src", "SpelledOut.jl")); using .SpelledOut
using Test
using SpelledOut

@testset "English" begin
    @test spelled_out(1234, lang = :en) == "one thousand, two hundred thirty-four"
    @test spelled_out(1234, lang = :en_US) == "one thousand, two hundred thirty-four"
    @test spelled_out(1234, lang = :en_UK) == "one thousand, two hundred and thirty-four"
    @test spelled_out(123456789, lang = :en) == "one hundred twenty-three million, four hundred fifty-six thousand, seven hundred eighty-nine"
    @test spelled_out(123456789, lang = :en_UK) == "one hundred and twenty-three million, four hundred and fifty-six thousand, seven hundred and eighty-nine"
    @test spelled_out(0, lang = :en) == "zero"
    @test spelled_out(0, lang = :en_UK) == "zero"
    @test spelled_out(100, lang = :en) == "one hundred"
    @test spelled_out(100, lang = :en_UK) == "one hundred"
    @test spelled_out(123456789, lang = :en_UK, dict=:european) == "one hundred and twenty-three million, four hundred and fifty-six thousand, seven hundred and eighty-nine"
    @test spelled_out(123456789, lang = :en, dict=:european) == "one hundred twenty-three million, four hundred fifty-six thousand, seven hundred eighty-nine"
    @test spelled_out(123456789, lang = :en_UK, dict=:british) == "one hundred and twenty-three million, four hundred and fifty-six thousand, seven hundred and eighty-nine"
    @test spelled_out(1234567890123, lang = :en_UK, dict=:british) == "one billion, two hundred and thirty-four thousand million, five hundred and sixty-seven million, eight hundred and ninety thousand, one hundred and twenty-three"
    @test spelled_out(1234567890123, lang = :en_UK, dict=:modern) == "one trillion, two hundred and thirty-four billion, five hundred and sixty-seven million, eight hundred and ninety thousand, one hundred and twenty-three"
    @test spelled_out(3.0, lang = :en) == "three"
    @test spelled_out(3, lang = :en) == "three"
    @test spelled_out(3.141592653589793, lang = :en) == "three point one four one five nine two six five three five eight nine seven nine three"
    @test spelled_out(2^log2(3390000), lang = :en) == "three million, three hundred ninety thousand" # should convert 3.3899999999999986e6 to an integer
    @test spelled_out(1000000.01, lang = :en) == "one million point zero one" # parse big floats correctly (i.e., avoid scientific notation.  Caveat: this is slow)
    @test spelled_out(1//3, lang = :en_UK) == "one third"
    @test spelled_out(1//5, lang = :en_UK) == "one fifth"
    @test spelled_out(1//102, lang = :en_UK) == "one one hundred and second"
    @test spelled_out(1//102, lang = :en_US) == "one one hundred second"
    @test spelled_out(5//102, lang = :en_UK) == "five one hundred and seconds"
    @test spelled_out(5//12, lang = :en_UK) == "five twelfths"
    @test spelled_out(1//1, lang = :en_UK) == "one"
    @test spelled_out(40//1, lang = :en_UK) == "forty"
    @test spelled_out(40//2, lang = :en_UK) == "twenty"
    @test spelled_out(40//6, lang = :en_UK) == "twenty thirds"
end

@testset "Spanish" begin
    @test spelled_out(585000, lang = :es) == "quinientos ochenta y cinco mil"

end

nothing