text:
default_ygrid(; lo=-6, hi=6, length=200) = range(lo, hi, length=length)

function group(sym::Symbol, chain)
  for monitor in chain
    if sym in keys(monitor[1])
      return [s[sym] for s in monitor]
    end
  end
end

function posterior_density(chain, ygrid)
  mu = group(:mu, chain)
  sigma = group(:sigma, chain)
  nu = group(:nu, chain)
  phi = group(:phi, chain)
  etaC = group(:etaC, chain)
  etaT = group(:etaT, chain)
  beta = group(:beta, chain)
  K = length(mu[1])
  B = length(beta)
  etaT_star = [beta[b] ? etaT[b] : etaC[b] for b in 1:B]

  function compute_pdf(b, eta)
    comps = [SkewT(mu[b][k], sigma[b][k], nu[b][k], phi[b][k]) for k in 1:K]
    pdf.(MixtureModel(comps , eta[b]), ygrid)
  end

  pdfC = [compute_pdf(b, etaC) for b in 1:B]
  pdfT = [compute_pdf(b, etaT_star) for b in 1:B]

  return (pdfC=pdfC, pdfT=pdfT)
end

function printsummary(chain, summarystats; laststate=nothing, digits=3)
  ll = [s[:loglike] for s in summarystats]
  p = group(:p, chain)
  gammaC = group(:gammaC, chain)
  gammaT = group(:gammaT, chain)
  phi = group(:phi, chain)
  etaC = group(:etaC, chain)
  etaT = group(:etaT, chain)
  mu = group(:mu, chain)
  beta = group(:beta, chain)
  sigma = group(:sigma, chain)
  phi = group(:phi, chain)
  nu = group(:nu, chain)
  psi = group(:psi, chain)
  omega = group(:omega, chain)
  tau = group(:tau, chain)

  sanitize(s) = round(mean(s), digits=digits)
  sanitizevec(v) = round.(mean(v), digits=digits)

  p == nothing || println("mean p: ", sanitize(p))
  println("mean beta: ", sanitize(beta))
  println("mean gammaC: ", sanitize(gammaC))
  println("mean gammaT: ", sanitize(gammaT))
  any(beta) && println("mean gammaT[β=1]: ", sanitize(gammaT[beta .== 1]))
  all(beta) || println("mean gammaT[β=0]: ", sanitize(gammaT[beta .== 0]))
  println("mean etaC: ", sanitizevec(etaC))
  println("mean etaT: ", sanitizevec(etaT))
  any(beta) && println("mean etaT[β=1]: ", sanitizevec(etaT[beta .== 1]))
  all(beta) || println("mean etaT[β=0]: ", sanitizevec(etaT[beta .== 0]))
  println("mean mu: ", sanitizevec(mu))
  println("mean sigma: ", sanitizevec(sigma))
  println("mean nu: ", sanitizevec(nu))
  println("mean phi: ", sanitizevec(phi))
  println("mean tau: ", sanitize(tau))

  if laststate != nothing
  end
end

function plot_posterior_predictive_only(chain; ygrid=default_ygrid(), lw=0,
                                        ls=:solid, labelCppd="C",
                                        labelTppd="T", alpha=0.3,
                                        density_legend_pos=:best,
                                        legendfont=font(12))
  pdfC, pdfT = posterior_density(chain, ygrid)

  pdfC_lower = MCMC.quantiles(hcat(pdfC...), .025, dims=2)
  pdfC_upper = MCMC.quantiles(hcat(pdfC...), .975, dims=2)
  pdfC_mean = mean(pdfC)

  pdfT_lower = MCMC.quantiles(hcat(pdfT...), .025, dims=2)
  pdfT_upper = MCMC.quantiles(hcat(pdfT...), .975, dims=2)
  pdfT_mean = mean(pdfT)

  # Plot density
  plot(ygrid, pdfC_lower, fillrange=pdfC_upper, alpha=alpha, color=:blue,
       label=labelCppd, legend=density_legend_pos)
  plot!(ygrid, pdfT_lower, fillrange=pdfT_upper, alpha=alpha, color=:red,
        label=labelTppd)
  if lw > 0
    plot!(ygrid, pdfC_mean, lw=lw, labels=nothing, color=:blue)
    plot!(ygrid, pdfT_mean, lw=lw, labels=nothing, color=:red)
  end
end

function plot_posterior_predictive(yC, yT, chain, bw; lw=.5, labelyC=L"\tilde y_C",
                                   ygrid=default_ygrid(), ls=:solid,
                                   labelyT=L"\tilde{y}_T", legendfontsize=12,
                                   alpha=0.3,
                                   labelCppd=nothing,
                                   labelTppd=nothing,
                                   simdata=nothing, density_legend_pos=:best)
  pdfC, pdfT = posterior_density(chain, ygrid)
  beta = group(:beta, chain)

  pdfC_lower = MCMC.quantiles(hcat(pdfC...), .025, dims=2)
  pdfC_upper = MCMC.quantiles(hcat(pdfC...), .975, dims=2)
  pdfC_mean = mean(pdfC)

  # if any(beta)
  pdfT_lower = MCMC.quantiles(hcat(pdfT...), .025, dims=2)
  pdfT_upper = MCMC.quantiles(hcat(pdfT...), .975, dims=2)
  pdfT_mean = mean(pdfT)
  # end

  # Plot density
  legendfont = font(12)

  if simdata == nothing
    plot(kde(yC[isfinite.(yC)], bandwidth=bw), lw=lw, label=labelyC,
         color=:blue, ls=ls) # legendfont=legendfont,
  else
    plot(ygrid, pdf.(simdata[:mmC], ygrid), lw=lw, label=labelyC,
         color=:blue, ls=ls)
  end
  plot!(ygrid, pdfC_lower, fillrange=pdfC_upper, alpha=alpha, color=:blue,
        label=labelCppd, legend=density_legend_pos)

  if simdata == nothing
    plot!(kde(yT[isfinite.(yT)], bandwidth=bw), lw=lw, label=labelyT,
          color=:red, ls=ls)
  else
    plot!(ygrid, pdf.(simdata[:mmT], ygrid), lw=lw, label=labelyT,
          color=:red, ls=ls)
  end
  # if any(beta)
  plot!(ygrid, pdfT_lower, fillrange=pdfT_upper, alpha=alpha, color=:red,
        label=labelTppd)
  # end
end

function plot_gamma(yC, yT, chain; hist=true, simdata=nothing, xrotation=nothing)
  gammaC = group(:gammaC, chain)
  gammaT = group(:gammaT, chain)
  beta = group(:beta, chain)

  # gammaT = any(beta) ? gammaT[beta] : gammaC
  gammaT_star = [beta[b] ? gammaT[b] : gammaC[b] for b in eachindex(beta)]

  if hist
    histogram(gammaC, color=:blue, normalize=true, alpha=0.3, linealpha=0,
              label=nothing, grid=false)
    histogram!(gammaT_star, color=:red, normalize=true, alpha=0.3, linealpha=0,
               label=nothing)
    xrotation == nothing || plot!(xrotation=xrotation)
    if simdata == nothing
      # Plot empirical mean
      vline!([mean(isinf.(yC)), mean(isinf.(yT))], color=[:blue, :red],
             label=nothing, ls=:dash)
    else
      vline!([simdata[:gammaC], simdata[:gammaT]], color=[:blue, :red],
             label=nothing, ls=:dash)
    end
    # Posterior mean.
    vline!([mean(gammaC), mean(gammaT_star)], color=[:blue, :red], alpha=.6,
           label=nothing)
    # 95% CI
    vline!([quantile(gammaC, .025), quantile(gammaC, .975)],
           color=:blue, label=nothing, ls=:dot, lw=2)
    vline!([quantile(gammaT_star, .025), quantile(gammaT_star, .975)],
           color=:red, label=nothing, ls=:dot, lw=2)

    xlabel!(L"\gamma_i")
    ylabel!("density")
  else
    boxplot(gammaC, outliers=false, color=:blue, label="", alpha=.5)
    boxplot!(gammaT_star, outliers=false, color=:red, label="", alpha=.5)
    xticks!([1, 2], [L"\gamma_C", L"\gamma_T"],
            xtickfont=font(20), ytickfont=font(16))
    scatter!([1, 2], [mean(isinf.(yC)), mean(isinf.(yT))], markersize=[10, 10],
             color=[:blue, :red], label=nothing)
  end
end

function boxplot_kernel_param(sym, chain; paramname="", simdata=nothing,
                              simsym=nothing)
  param = group(sym, chain)
  boxplot(hcat(param...)', outliers=false, labels=nothing)
  xlabel!("$(paramname) Mixture Components", font=font(12))
  simdata == nothing || hline!(simdata[simsym], ls=:dot, lw=3, labels="truth")
  return
end

function trace_kernel_mean(sym, chain; paramname="", simdata=nothing,
                           simsym=nothing)
  param = group(sym, chain)
  nsamps = length(param)
  plot(cumsum(param) ./ collect(1:nsamps), outliers=false, labels=nothing)
  xlabel!("MCMC iteration", font=font(12))
  ylabel!("Posterior Expectation of $(paramname)", font=font(12))
  simdata == nothing || hline!(simdata[simsym], ls=:dot, labels="truth")
  return
end

function trace_kernel_param(sym, chain; paramname="", simdata=nothing,
                            simsym=nothing)
  param = group(sym, chain)
  plot(hcat(param...)', outliers=false, labels=nothing)
  xlabel!("MCMC iteration", font=font(12))
  ylabel!(paramname, font=font(12))
  simdata == nothing || hline!(simdata[simsym], ls=:dot, labels="truth")
  return
end

function plot_gamma(yC, yT, chain0, chain1, pm1, imgdir; plotsize=(400, 400),
                    hist=true, simdata=nothing, xrotation=nothing)
  @assert length(chain0[1]) == length(chain1[1])
  B = length(chain0[1])
  chain = [[pm1 > rand() ? chain1[1][b] : chain0[1][b] for b in 1:B]]

  plot_gamma(yC, yT, chain, hist=hist, simdata=simdata)
  plot!(size=plotsize)
  xrotation == nothing || plot!(xrotation=xrotation)
  savefig("$(imgdir)/gammas.pdf")
  closeall()
end

function plot_posterior_predictive(yC, yT, chain0, chain1, pm1, imgdir;
                                   bw_postpred=0.3, simdata=nothing,
                                   ygrid=default_ygrid(), xlims_=nothing,
                                   plotsize=(400,400), lw=.5, ls=:solid,
                                   density_legend_pos=:best, showpm1=false,
                                   digits=3, fontsize=12,
                                   binsC=:auto, binsT=:auto)
  @assert length(chain0[1]) == length(chain1[1])
  B = length(chain0[1])
  chain = [[pm1 > rand() ? chain1[1][b] : chain0[1][b] for b in 1:B]]

  pm1 = round(pm1, digits=digits)

  # Posterior predictive density
  plot_posterior_predictive(yC, yT, chain, bw_postpred, ygrid=ygrid, ls=ls,
                            density_legend_pos=density_legend_pos, lw=lw)
  xlims_ == nothing || xlims!(xlims_)
  plot!(size=plotsize)
  showpm1 && plot!(legendtitle="P(M₁|y) ≈ $pm1", legendtitlefont=font(fontsize))
  savefig("$(imgdir)/postpred.pdf")
  closeall()

  # Posterior predictive density and data histogram
  plot_posterior_predictive_only(chain, ygrid=ygrid, lw=.5,
                                 density_legend_pos=density_legend_pos)
  plot!(size=plotsize)
  histogram!(yC, alpha=0.3, linealpha=0, normalize=true, labels=nothing,
             color=:grey, bins=binsC)
  histogram!(yT, alpha=0.3, linealpha=0, normalize=true, labels=nothing,
             color=:grey, bins=binsT)
  savefig("$(imgdir)/postpred-data-hist.pdf")
  closeall()

  # Posterior predictive density, mean, and sim-truth.
  # TODO: check
  if simdata != nothing
    plot_posterior_predictive_only(chain, ygrid=ygrid, lw=.5,
                                   density_legend_pos=density_legend_pos)
    plot!(ygrid, pdf.(simdata[:mmC], ygrid), lw=1, label=nothing,
         color=:blue, ls=:dot)
    plot!(ygrid, pdf.(simdata[:mmT], ygrid), lw=1, label=nothing,
          color=:red, ls=:dot)
    plot!(size=plotsize)
    savefig("$(imgdir)/postpred-mean-simtruth.pdf")
    closeall()
  end

  # Posterior predictive density with true data density
  if simdata != nothing
    plot_posterior_predictive(yC, yT, chain, bw_postpred, ygrid=ygrid,
                              simdata=simdata, lw=lw, ls=ls,
                              density_legend_pos=density_legend_pos)
    xlims_ == nothing || xlims!(xlims_)
    plot!(size=plotsize)
    showpm1 && plot!(legendtitle="P(M₁|y) ≈ $pm1", legendtitlefont=font(fontsize))
    savefig("$(imgdir)/postpred-true-data-density.pdf")
    closeall()
  end
end

function plotpostsummary(chain, summarystats, yC, yT, imgdir; digits=3,
                         laststate=nothing, bw_postpred=0.3, simdata=nothing,
                         ygrid=default_ygrid(), xlims_=nothing, gamma_hist=true,
                         plotsize=(400,400), density_legend_pos=:best)
  # Loglikelihood trace
  ll = [s[:loglike] for s in summarystats]
  plot(ll, label=nothing)
  ylabel!("log likelihood")
  xlabel!("MCMC iteration")
  plot!(size=plotsize)
  savefig("$(imgdir)/loglike.pdf")
  closeall()

  # Posterior density
  plot_posterior_predictive(yC, yT, chain, bw_postpred, ygrid=ygrid,
                            density_legend_pos=density_legend_pos)
  xlims_ == nothing || xlims!(xlims_)
  plot!(size=plotsize)
  savefig("$(imgdir)/postpred.pdf")
  closeall()

  if simdata != nothing
    plot_posterior_predictive(yC, yT, chain, bw_postpred, ygrid=ygrid,
                              simdata=simdata,
                              density_legend_pos=density_legend_pos)
    xlims_ == nothing || xlims!(xlims_)
    plot!(size=plotsize)
    savefig("$(imgdir)/postpred-true-data-density.pdf")
    closeall()
  end

  if haskey(chain[1][1], :tau)
    plot(group(:tau, chain), label=nothing, size=plotsize)
    xlabel!("MCMC iteration", font=font(12))
    ylabel!("τ", font=font(12))
    savefig("$(imgdir)/tau-trace.pdf"); closeall()
  end

  # Trace of p, beta.
  if haskey(chain[1][1], :p)
    trace_kernel_param(:p, chain, paramname="p")
    plot!(size=plotsize)
    savefig("$(imgdir)/p-trace.pdf"); closeall()

    trace_kernel_mean(:p, chain, paramname="p")
    plot!(size=plotsize)
    savefig("$(imgdir)/p-mean-trace.pdf"); closeall()
  end
  trace_kernel_param(:beta, chain, paramname="β")
  plot!(size=plotsize)
  savefig("$(imgdir)/beta-trace.pdf"); closeall()
  trace_kernel_mean(:beta, chain, paramname="β")
  plot!(size=plotsize)
  hline!([0.5], ls=:dot, labels=nothing)
  savefig("$(imgdir)/beta-mean-trace.pdf"); closeall()

  # Proportion of gammas.
  plot_gamma(yC, yT, chain, hist=gamma_hist)
  plot!(size=plotsize)
  savefig("$(imgdir)/gammas.pdf")
  closeall()

  # Plot kernel parameters boxplot
  boxplot_kernel_param(:mu, chain, paramname="μ", simdata=simdata, simsym=:loc)
  plot!(size=plotsize)
  savefig("$(imgdir)/mu.pdf"); closeall()

  boxplot_kernel_param(:sigma, chain, paramname="σ", simdata=simdata, simsym=:scale)
  plot!(size=plotsize)
  savefig("$(imgdir)/sigma.pdf"); closeall()

  boxplot_kernel_param(:nu, chain, paramname="ν", simdata=simdata, simsym=:df)
  plot!(size=plotsize)
  savefig("$(imgdir)/nu.pdf"); closeall()

  boxplot_kernel_param(:phi, chain, paramname="ϕ", simdata=simdata, simsym=:skew)
  plot!(size=plotsize)
  savefig("$(imgdir)/phi.pdf"); closeall()

  boxplot_kernel_param(:psi, chain, paramname="ψ")
  plot!(size=plotsize)
  savefig("$(imgdir)/psi.pdf"); closeall()

  boxplot_kernel_param(:omega, chain, paramname="ω")
  plot!(size=plotsize)
  savefig("$(imgdir)/omega.pdf"); closeall()

  boxplot_kernel_param(:etaC, chain, paramname="ηc", simdata=simdata, simsym=:etaC)
  plot!(size=plotsize)
  savefig("$(imgdir)/etaC.pdf"); closeall()

  boxplot_kernel_param(:etaT, chain, paramname="ηt", simdata=simdata, simsym=:etaT)
  plot!(size=plotsize)
  savefig("$(imgdir)/etaT.pdf"); closeall()

  # Plot kernel parameters trace plot
  trace_kernel_param(:mu, chain, paramname="μ")
  plot!(size=plotsize)
  savefig("$(imgdir)/mu-trace.pdf"); closeall()

  trace_kernel_param(:sigma, chain, paramname="σ")
  plot!(size=plotsize)
  savefig("$(imgdir)/sigma-trace.pdf"); closeall()

  trace_kernel_param(:nu, chain, paramname="ν")
  plot!(size=plotsize)
  savefig("$(imgdir)/nu-trace.pdf"); closeall()

  trace_kernel_param(:phi, chain, paramname="ϕ")
  plot!(size=plotsize)
  savefig("$(imgdir)/phi-trace.pdf"); closeall()

  trace_kernel_param(:etaC, chain, paramname="ηc")
  plot!(size=plotsize)
  savefig("$(imgdir)/etaC-trace.pdf"); closeall()

  trace_kernel_param(:etaT, chain, paramname="ηt")
  plot!(size=plotsize)
  savefig("$(imgdir)/etaT-trace.pdf"); closeall()

  trace_kernel_param(:psi, chain, paramname="ψ")
  plot!(size=plotsize)
  savefig("$(imgdir)/psi-trace.pdf"); closeall()

  trace_kernel_param(:omega, chain, paramname="ω")
  plot!(size=plotsize)
  savefig("$(imgdir)/omega-trace.pdf"); closeall()
end

max_stars_repo_name: luiarthur/CytofDensityEstimation

max_issues_repo_name: luiarthur/CytofDensityEstimation

max_forks_repo_name: luiarthur/CytofDensityEstimation


text:
module RocketTest

using Test, Documenter, Rocket

doctest(Rocket)

@testset "Rocket" begin

    include("./test_teardown.jl")
    include("./teardown/test_void_teardown.jl")

    include("./test_actor.jl")
    include("./actor/test_test_actor.jl")
    include("./actor/test_void_actor.jl")
    include("./actor/test_lambda_actor.jl")
    include("./actor/test_logger_actor.jl")
    include("./actor/test_keep_actor.jl")
    include("./actor/test_buffer_actor.jl")
    include("./actor/test_sync_actor.jl")
    include("./actor/test_function_actor.jl")

    include("./test_subscribable.jl")
    include("./observable/test_observable_function.jl")
    include("./observable/test_observable_single.jl")
    include("./observable/test_observable_array.jl")
    include("./observable/test_observable_iterable.jl")
    include("./observable/test_observable_faulted.jl")
    include("./observable/test_observable_never.jl")
    include("./observable/test_observable_completed.jl")
    include("./observable/test_observable_timer.jl")
    include("./observable/test_observable_interval.jl")
    include("./observable/test_observable_proxy.jl")
    include("./observable/test_observable_combine_latest.jl")
    include("./observable/test_observable_combine_updates.jl")
    include("./observable/test_observable_collect_latest.jl")
    include("./observable/test_observable_race.jl")
    include("./observable/test_observable_merged.jl")
    include("./observable/test_observable_concat.jl")
    include("./observable/test_observable_generate.jl")
    include("./observable/test_observable_defer.jl")
    include("./observable/test_observable_zipped.jl")

    include("./test_subject.jl")
    include("./subjects/test_subject.jl")
    include("./subjects/test_subject_behavior.jl")
    include("./subjects/test_subject_replay.jl")
    include("./subjects/test_subject_pending.jl")

    include("./test_operator.jl")
    include("./operators/test_operator_async.jl")
    include("./operators/test_operator_map.jl")
    include("./operators/test_operator_map_to.jl")
    include("./operators/test_operator_scan.jl")
    include("./operators/test_operator_enumerate.jl")
    include("./operators/test_operator_uppercase.jl")
    include("./operators/test_operator_lowercase.jl")
    include("./operators/test_operator_to_array.jl")
    include("./operators/test_operator_filter.jl")
    include("./operators/test_operator_find.jl")
    include("./operators/test_operator_find_index.jl")
    include("./operators/test_operator_some.jl")
    include("./operators/test_operator_take.jl")
    include("./operators/test_operator_take_until.jl")
    include("./operators/test_operator_first.jl")
    include("./operators/test_operator_last.jl")
    include("./operators/test_operator_count.jl")
    include("./operators/test_operator_max.jl")
    include("./operators/test_operator_min.jl")
    include("./operators/test_operator_reduce.jl")
    include("./operators/test_operator_sum.jl")
    include("./operators/test_operator_catch_error.jl")
    include("./operators/test_operator_rerun.jl")
    include("./operators/test_operator_tap.jl")
    include("./operators/test_operator_tap_on_subscribe.jl")
    include("./operators/test_operator_tap_on_complete.jl")
    include("./operators/test_operator_delay.jl")
    include("./operators/test_operator_safe.jl")
    include("./operators/test_operator_noop.jl")
    include("./operators/test_operator_multicast.jl")
    include("./operators/test_operator_ref_count.jl")
    include("./operators/test_operator_default_if_empty.jl")
    include("./operators/test_operator_error_if_empty.jl")
    include("./operators/test_operator_skip_next.jl")
    include("./operators/test_operator_skip_error.jl")
    include("./operators/test_operator_skip_complete.jl")
    include("./operators/test_operator_tuple_with.jl")
    include("./operators/test_operator_switch_map.jl")
    include("./operators/test_operator_switch_map_to.jl")
    include("./operators/test_operator_merge_map.jl")
    include("./operators/test_operator_concat_map.jl")
    include("./operators/test_operator_concat_map_to.jl")
    include("./operators/test_operator_exhaust_map.jl")
    include("./operators/test_operator_discontinue.jl")
    include("./operators/test_operator_with_latest.jl")
    include("./operators/test_operator_ignore.jl")
    include("./operators/test_operator_start_with.jl")
    include("./operators/test_operator_accumulated.jl")
    include("./operators/test_operator_pairwise.jl")

    @testset "Detect ambiguities" begin
        @test length(Test.detect_ambiguities(Rocket)) == 0
    end
end

end

ext: jl


text:
using NonconvexBayesian, NonconvexIpopt, NonconvexNLopt
using LinearAlgebra, Test

f(x::AbstractVector) = sqrt(x[2])
g(x::AbstractVector, a, b) = (a*x[1] + b)^3 - x[2]

@testset "Cheap objective and constraints" begin
    @testset "Fit prior: $fit_prior" for fit_prior in (true, false)
        m = Model()
        set_objective!(m, f)
        addvar!(m, [1e-4, 1e-4], [10.0, 10.0])
        add_ineq_constraint!(m, x -> g(x, 2, 0))
        add_ineq_constraint!(m, x -> g(x, -1, 1))

        alg = BayesOptAlg(IpoptAlg())
        options = BayesOptOptions(
            sub_options = IpoptOptions(print_level = 0), maxiter = 10, ctol = 1e-4,
            ninit = 2, initialize = true, postoptimize = false, fit_prior = fit_prior,
        )
        r = optimize(m, alg, [1.234, 2.345], options = options)
        @test abs(r.minimum - sqrt(8/27)) < 1e-3
        @test norm(r.minimizer - [1/3, 8/27]) < 1e-3
    end
end

@testset "Expensive objective" begin
    @testset "Fit prior: $fit_prior" for fit_prior in (true, false)
        m = Model()
        set_objective!(m, f, flags = [:expensive])
        addvar!(m, [1e-4, 1e-4], [10.0, 10.0])
        add_ineq_constraint!(m, x -> g(x, 2, 0))
        add_ineq_constraint!(m, x -> g(x, -1, 1))
        add_eq_constraint!(m, x -> sum(x) - 1/3 - 8/27)

        alg = BayesOptAlg(IpoptAlg())
        options = BayesOptOptions(
            sub_options = IpoptOptions(print_level = 0), maxiter = 10, ctol = 1e-4,
            ninit = 2, initialize = true, postoptimize = false, skip = 2, every = 2,
            fit_prior = fit_prior,
        )
        r = optimize(m, alg, [1.234, 2.345], options = options)
        @test abs(r.minimum - sqrt(8/27)) < 1e-3
        @test norm(r.minimizer - [1/3, 8/27]) < 1e-3
    end
end

@testset "Expensive objective - postoptimize" begin
    @testset "Fit prior: $fit_prior" for fit_prior in (true, false)
        m = Model()
        set_objective!(m, f, flags = [:expensive])
        addvar!(m, [1e-4, 1e-4], [10.0, 10.0])
        add_ineq_constraint!(m, x -> g(x, 2, 0))
        add_ineq_constraint!(m, x -> g(x, -1, 1))
        add_eq_constraint!(m, x -> sum(x) - 1/3 - 8/27)

        alg = BayesOptAlg(IpoptAlg())
        options = BayesOptOptions(
            sub_options = IpoptOptions(print_level = 0), maxiter = 10, ctol = 1e-4,
            ninit = 2, initialize = true, postoptimize = true, fit_prior = fit_prior,
        )
        r = optimize(m, alg, [1.234, 2.345], options = options)
        @test abs(r.minimum - sqrt(8/27)) < 1e-5
        @test norm(r.minimizer - [1/3, 8/27]) < 1e-5
    end
end

@testset "Expensive objective - no preinitialize" begin
    @testset "Fit prior: $fit_prior" for fit_prior in (true, false)
        m = Model()
        set_objective!(m, f, flags = [:expensive])
        addvar!(m, [1e-4, 1e-4], [10.0, 10.0])
        add_ineq_constraint!(m, x -> g(x, 2, 0))
        add_ineq_constraint!(m, x -> g(x, -1, 1))
        add_eq_constraint!(m, x -> sum(x) - 1/3 - 8/27)

        alg = BayesOptAlg(IpoptAlg())
        options = BayesOptOptions(
            sub_options = IpoptOptions(print_level = 0), maxiter = 10, ctol = 1e-4,
            initialize = false, postoptimize = false, skip = 2, every = 2, fit_prior = fit_prior,
        )
        r = optimize(m, alg, [1.234, 2.345], options = options)
        @test abs(r.minimum - sqrt(8/27)) < 1e-5
        @test norm(r.minimizer - [1/3, 8/27]) < 1e-5
    end
end

@testset "Expensive inequality constraint" begin
    @testset "Fit prior: $fit_prior" for fit_prior in (true, false)
        m = Model()
        set_objective!(m, f)
        addvar!(m, [1e-4, 1e-4], [10.0, 10.0])
        add_ineq_constraint!(m, x -> g(x, 2, 0), flags = [:expensive])
        add_ineq_constraint!(m, x -> g(x, -1, 1))

        alg = BayesOptAlg(IpoptAlg())
        options = BayesOptOptions(
            sub_options = IpoptOptions(print_level = 0), maxiter = 100, ctol = 1e-4,
            ninit = 5, initialize = true, postoptimize = false, last = 20, fit_prior = fit_prior,
        )
        r = optimize(m, alg, [1.234, 2.345], options = options)
        @test abs(r.minimum - sqrt(8/27)) < 1e-3
        @test norm(r.minimizer - [1/3, 8/27]) < 1e-3
    end
end

@testset "Expensive inequality constraint with equality" begin
    @testset "Fit prior: $fit_prior" for fit_prior in (true, false)
        m = Model()
        set_objective!(m, f)
        addvar!(m, [1e-4, 1e-4], [10.0, 10.0])
        add_ineq_constraint!(m, x -> g(x, 2, 0), flags = [:expensive])
        add_ineq_constraint!(m, x -> g(x, -1, 1))

        alg = BayesOptAlg(IpoptAlg())
        options = BayesOptOptions(
            sub_options = IpoptOptions(print_level = 0), maxiter = 20, ctol = 1e-4,
            ninit = 2, initialize = true, postoptimize = false, last = 10, fit_prior = fit_prior,
        )
        r = optimize(m, alg, [1.234, 2.345], options = options)
        @test abs(r.minimum - sqrt(8/27)) < 1e-3
        @test norm(r.minimizer - [1/3, 8/27]) < 1e-3
    end
end

@testset "Two expensive inequality constraint" begin
    @testset "Fit prior: $fit_prior" for fit_prior in (true, false)
        m = Model()
        set_objective!(m, f)
        addvar!(m, [1e-4, 1e-4], [10.0, 10.0])
        add_ineq_constraint!(m, x -> g(x, 2, 0), flags = [:expensive])
        add_ineq_constraint!(m, x -> g(x, -1, 1), flags = [:expensive])

        alg = BayesOptAlg(IpoptAlg())
        options = BayesOptOptions(
            sub_options = IpoptOptions(print_level = 0), maxiter = 40, ctol = 1e-4,
            ninit = 2, initialize = true, postoptimize = false, every = 5, skip = 2,
            last = 10, fit_prior = fit_prior,
        )
        r = optimize(m, alg, [1.234, 2.345], options = options)
        @test abs(r.minimum - sqrt(8/27)) < 1e-3
        @test norm(r.minimizer - [1/3, 8/27]) < 1e-3
    end
end

@testset "Expensive equality constraint" begin
    @testset "Fit prior: $fit_prior" for fit_prior in (true, false)
        m = Model()
        set_objective!(m, f)
        addvar!(m, [1e-4, 1e-4], [10.0, 10.0])
        add_ineq_constraint!(m, x -> g(x, 2, 0))
        add_ineq_constraint!(m, x -> g(x, -1, 1))
        add_eq_constraint!(m, x -> sum(x) - 1/3 - 8/27, flags = [:expensive])
        alg = BayesOptAlg(IpoptAlg())
        options = BayesOptOptions(
            sub_options = IpoptOptions(print_level = 0), maxiter = 30, ctol = 1e-4,
            ninit = 2, initialize = true, postoptimize = false, last = 10, fit_prior = fit_prior,
        )
        r = optimize(
            m, alg, [1.234, 2.345], options = options,
        )
        @test abs(r.minimum - sqrt(8/27)) < 1e-3
        @test norm(r.minimizer - [1/3, 8/27]) < 1e-3
    end
end

@testset "Expensive block constraint" begin
    @testset "Fit prior: $fit_prior" for fit_prior in (true, false)
        m = Model()
        set_objective!(m, f)
        addvar!(m, [1e-4, 1e-4], [10.0, 10.0])
        add_ineq_constraint!(m, x -> [g(x, 2, 0), g(x, -1, 1)], flags = [:expensive])
        alg = BayesOptAlg(IpoptAlg())
        options = BayesOptOptions(
            sub_options = IpoptOptions(print_level = 0), maxiter = 20, ctol = 1e-4,
            ninit = 2, initialize = true, postoptimize = false, skip = 2, every = 2,
            last = 10, fit_prior = fit_prior,
        )
        r = optimize(
            m, alg, [1.234, 2.345], options = options,
        )
        @test abs(r.minimum - sqrt(8/27)) < 1e-3
        @test norm(r.minimizer - [1/3, 8/27]) < 1e-3
    end
end

@testset "All expensive no equality" begin
    @testset "Fit prior: $fit_prior" for fit_prior in (true, false)
        m = Model()
        set_objective!(m, f, flags = [:expensive])
        addvar!(m, [1e-4, 1e-4], [10.0, 10.0])
        add_ineq_constraint!(m, x -> [g(x, 2, 0), g(x, -1, 1)], flags = [:expensive])
        alg = BayesOptAlg(IpoptAlg())
        options = BayesOptOptions(
            sub_options = IpoptOptions(), maxiter = 50, ctol = 1e-3,
            ninit = 2, initialize = true, postoptimize = false, last = 10,
            fit_prior = fit_prior,
        )
        r = optimize(m, alg, [1.234, 2.345], options = options)
        @test abs(r.minimum - sqrt(8/27)) < 1e-3
        @test norm(r.minimizer - [1/3, 8/27]) < 1e-3
    end
end

@testset "All expensive with equality" begin
    @testset "Fit prior: $fit_prior" for fit_prior in (true, false)
        m = Model()
        set_objective!(m, f, flags = [:expensive])
        addvar!(m, [1e-4, 1e-4], [10.0, 10.0])
        add_ineq_constraint!(m, x -> [g(x, 2, 0), g(x, -1, 1)], flags = [:expensive])
        add_eq_constraint!(m, x -> sum(x) - 1/3 - 8/27, flags = [:expensive])
        alg = BayesOptAlg(IpoptAlg())
        options = BayesOptOptions(
            sub_options = IpoptOptions(print_level = 0), maxiter = 100, ctol = 1e-4,
            ninit = 2, initialize = true, postoptimize = false, last = 30, skip = 2,
            every = 2, fit_prior = fit_prior,
        )
        r = optimize(m, alg, [1.234, 2.345], options = options)
        @test abs(r.minimum - sqrt(8/27)) < 1e-3
        @test norm(r.minimizer - [1/3, 8/27]) < 1e-3
    end
end

@testset "Passing surrogates in" begin
    @testset "Fit prior: $fit_prior" for fit_prior in (true, false)
        x0 = [1.234, 2.345]
        s1 = NonconvexBayesian.surrogate(f, x0)
        s2 = NonconvexBayesian.surrogate(x -> [g(x, 2, 0), g(x, -1, 1)], x0)

        m = Model()
        set_objective!(m, x -> s1(x).lo, flags = [:expensive])
        addvar!(m, [1e-4, 1e-4], [10.0, 10.0])
        add_ineq_constraint!(
            m, x -> getproperty.(s2(x), :lo), flags = [:expensive],
        )
        alg = BayesOptAlg(IpoptAlg())
        options = BayesOptOptions(
            sub_options = IpoptOptions(print_level = 0), maxiter = 50, ctol = 1e-4,
            ninit = 2, initialize = true, postoptimize = false, fit_prior = fit_prior,
        )
        r = optimize(m, alg, [1.234, 2.345], options = options, surrogates = [s1, s2])
        @test abs(r.minimum - sqrt(8/27)) < 1e-3
        @test norm(r.minimizer - [1/3, 8/27]) < 1e-3
    end
end


text:
# Blockchain.jl copyright © 2020 Derik Kauffman
# Based on https://hackernoon.com/learn-blockchains-by-building-one-117428612f46
# Original code copyright © 2017 Daniel van Flymen
# Licensed under the MIT license.

module Snowcoin

using SHA: sha256

struct Transaction
    sender::String
    recipient::String
    amount::Int
end

struct Block
    index::Int
    timestamp::Int  # The rough number of seconds since the Unix epoch
    transactions::Vector{Transaction}
    proof::Int
    previous_hash::String
end

mutable struct Blockchain
    chain::Vector{Block}
    pending_transactions::Vector{Transaction}  # "Mempool"
end

"Create a new blockchain with a genesis block"
function Blockchain()
    bc = Blockchain([], [])
	new_block!(bc, 0, shash(rand(Int)))
	return bc
end

Base.length(bc::Blockchain) = length(bc.chain)

function is_valid(bc::Blockchain)
	prev_block = bc.chain[1]
	curr_index = 2

	while curr_index <= length(bc)
		block = bc.chain[curr_index]
		prev_block_hash = shash(prev_block)

		# Check that the hash of the block is correct
		if block.previous_hash != prev_block_hash
			return false
		end

		# Check that the proof of work is correct
		if !valid_proof(prev_block.proof, block.proof, prev_block_hash)
			return false
		end

		prev_block = block
		curr_index += 1
	end

	return true
end

"Create a new block and add it to the chain."
function new_block!(bc::Blockchain, proof, previous_hash)
    block = Block(
        length(bc) + 1,
        round(time()),
        bc.pending_transactions,
        proof,
        previous_hash,
    )

    bc.pending_transactions = []

    push!(bc.chain, block)
    return block
end
new_block!(bc::Blockchain, proof) = new_block!(bc, proof, shash(bc.chain[end]))

"Create a new transaction to go into the next mined block."
function new_transaction!(bc::Blockchain, sender, recepient, amount)
    transaction = Transaction(sender, recepient, amount)
    push!(bc.pending_transactions, transaction)

    return bc.chain[end].index + 1
end

# Double hash to prevent against length extension
# TODO: Consider using a different (memory-hard) hash algorithm
#SHA.sha256(x::Integer) = sha256(reinterpret(UInt8, [x]))
shash(x) = x |> string |> sha256 |> sha256 |> bytes2hex

# The number of zero bits that a hash must begin with to constitute a valid proof.
# TODO: Adjust dynamically to hold block time approximately constant.
const DIFFICULTY = UInt8(12)

# The number of snowcoins generated upon successfully mining a block.
# TODO: Adjust dynamically to limit total currency supply.
const REWARD = 120

# The number of units (flakes) each snowcoin can be subdivided into.
const PRECISION = 100_000_000

function valid_proof(last_proof, proof, last_hash)
    guess = string(last_proof) * string(proof) * last_hash
    guess_hash = shash(guess)
    return guess_hash <= -1 >>> DIFFICULTY  # True if the first `difficulty` bits are 0
end

function proof_of_work(last_proof, last_hash)
    # The address in the generation transaction provides enough randomness
	# for the hash, so we can start checking proofs at zero.
	proof = 0

    while !valid_proof(last_proof, proof, last_hash)
        proof += 1
    end
    return proof
end

include("server.jl")

new_transaction!(snowchain, "Derik", "Andrew", 10)
print(snowchain)

end

max_stars_repo_name: derikk/Snowcoin

max_issues_repo_name: derikk/Snowcoin

max_forks_repo_name: derikk/Snowcoin


text:
# Recreate plot with normalized density vs. orbital period:

include("../../../src/regress.jl")
using Statistics
using PyPlot
using DelimitedFiles

period = [1.510826, 2.421937,4.049219, 6.101013, 9.207540,12.352446,18.772866]

dir = "../../../data/POSTERIOR_NORM_DENSITY/"
fname = ["out_T1b_norm","out_T1c_norm","out_T1d_norm","out_T1e_norm","out_T1f_norm","out_T1g_norm","out_T1h_norm"]

np = 7; nsamp = 10000

dens_norm = zeros(np,nsamp)
dens_mean = zeros(np)
dens_sig = zeros(np)
for ip=1:np
  data = readdlm(string(dir,fname[ip]))
  dens_norm[ip,:] .= data[:,2]
  dens_mean[ip] = mean(dens_norm[ip,:])
  dens_sig[ip] = std(dens_norm[ip,:])
end

# Now, carry out regressions for each of these:
coeff_samp = zeros(2,nsamp)
fn = zeros(2,np)
fn[1,:] .= 1.0
fn[2,:] .= period
sig = ones(np)
dens_avg = 0.0
for isamp=1:nsamp
  global dens_avg += mean(dens_norm[:,isamp])
  coeff,cov = regress(fn,dens_norm[:,isamp],sig)
  coeff_samp[:,isamp] = coeff
end
dens_avg /= nsamp
# Now, compute standard deviation:
dens_var = 0.0
for isamp=1:nsamp
  global dens_var += (mean(dens_norm[:,isamp]) - dens_avg)^2
end
dens_sig = sqrt(dens_var/nsamp)

# Now, draw 1-sigma and 2-sigma confidence intervals:
period_trial = collect(0.0:0.1:20.0)
conf = zeros(5,length(period_trial))
for (i,p) in enumerate(period_trial)
  # Predict samples at this period
  csamp = coeff_samp[1,:] .+ coeff_samp[2,:] .* p
  # Now, sort these and pick out 68.3% & 95% confidence limits:
  csamp = sort(csamp)
  conf[:,i] = [csamp[250],csamp[1415],csamp[5000],csamp[8415],csamp[9750]]
end
# Now, plot median, 1-sigma and 2-sigma confidence intervals:

clf()
#plot(period_trial,conf[1,:],color="C1",alpha=0.5,linestyle="--")
errorbar(period,dens_mean,yerr=dens_sig,fmt="o")
#plot(period_trial,conf[5,:],color="C1",alpha=0.5,linestyle="--")
plot(period_trial,conf[3,:],color="C1",linewidth=2)
plot(period_trial,conf[2,:],color="C1",linewidth=2,linestyle=":")
plot(period_trial,conf[4,:],color="C1",linewidth=2,linestyle=":")
#plot(period_trial,dens_avg  .+ 0.0 .* period_trial,color="C0",linewidth=2)
fill_between([0,20],[dens_avg-dens_sig,dens_avg-dens_sig],[dens_avg+dens_sig,dens_avg+dens_sig],alpha=0.2,color="C0")
axis([0,20,0,1.5])
xlabel("Orbital Period [d]")
ylabel("Normalized density")
tight_layout()
savefig("../Norm_dens_vs_period.pdf",bbox_inches = "tight")


text:
# =================================================================
# Van der Pol model
# See https://easychair.org/publications/paper/gjfh
# =================================================================

using Reachability, MathematicalSystems, TaylorIntegration

# Equations of motion: we write the function such that the operations are either
# unary or binary
@taylorize function vanderPol_mu_one!(dx, x, params, t)
    local μ = 1.0
    dx[1] = x[2]
    dx[2] = (μ * x[2]) * (1 - x[1]^2) - x[1]
    return dx
end

@taylorize function vanderPol_mu_two!(dx, x, params, t)
    local μ = 2.0
    dx[1] = x[2]
    dx[2] = (μ * x[2]) * (1 - x[1]^2) - x[1]
    return dx
end

function vanderpol(; μ=1.0,
                     X0=Hyperrectangle(low=[1.25, 2.35], high=[1.55, 2.45]),
                     T=7.0,
                     property=(t, x) -> x[2] < 2.75)

    if μ == 1.0
        f = vanderPol_mu_one!
    elseif μ == 2.0
        f = vanderPol_mu_two!
    else
        error("the value of μ = $μ is not implemented")
    end
    # equations, x' = f(x(t))
    F = BlackBoxContinuousSystem(f, 2)

    # problem options
    𝑂 = Options(:T=>T, :mode=>"check", :property=>property)

    # instantiate problem
    𝑃 = InitialValueProblem(F, X0)

    return 𝑃, 𝑂
end

function split(H::AbstractHyperrectangle, n::Int=2, m::Int=2)
     @assert dim(H) == 2
     r = copy(radius_hyperrectangle(H))
     r[1] /= n
     r[2] /= m
     d = 2*r
     result = Vector{Hyperrectangle}()
     sizehint!(result, n*m)
     c0 = low(H) - r
     c = copy(c0)
     for i in 1:n
         c[1] += d[1]
         c[2] = c0[2]
         for j in 1:m
             c[2] += d[2]
             push!(result, Hyperrectangle(copy(c), copy(r)))
         end
     end
     return result
end

max_stars_repo_name: JuliaReach/ReachabilityBenchmarks

max_issues_repo_name: JuliaReach/ReachabilityBenchmarks

max_forks_repo_name: JuliaReach/ReachabilityBenchmarks


text:
@testset "AAE" begin
	using DrWatson
	@quickactivate
	using Test
	using GenerativeAD
	import StatsBase: fit!, predict
	using ConditionalDists, GenerativeModels, Distributions
	using EvalMetrics
	using IPMeasures
	using Flux

	# toy example
	M = 4
	N1 = 100
	N2 = 50
	X = vcat(randn(Float32, Int(M/2), N1), ones(Float32, Int(M/2), N1))
	Y = vcat(randn(Float32, Int(M/2), N2), ones(Float32, Int(M/2), N2))
	Z = vcat(ones(Float32, Int(M/2), N2), randn(Float32, Int(M/2), N2))
	data = ((X,),(Y,zeros(size(Y,2))))
	parameters = (zdim=8, hdim=32, lr=0.001, batchsize=8, activation="swish", nlayers=3, idim=M)
	function basic_convergence_test(model)
		history, iterations, model = fit!(model, data; patience = 100, max_train_time=600, parameters...)
		scores = map(x->GenerativeAD.Models.reconstruction_score(model, x), (X,Y,Z))
		@test mean(scores[1]) < mean(scores[3])
		@test mean(scores[2]) < mean(scores[3])
		model
	end
	# AAE
	model = GenerativeAD.Models.aae_constructor(;parameters...)
	@test length(Flux.params(model)) == length(Flux.params(model.encoder)) +
		length(Flux.params(model.decoder)) + length(Flux.params(model.discriminator))
	model = basic_convergence_test(model)
	@test size(GenerativeAD.Models.reconstruct(model,X)) == size(X)
	@test length(GenerativeAD.Models.reconstruction_score_mean(model,X)) == size(X,2)
	@test length(GenerativeAD.Models.latent_score(model,X)) == size(X,2)
	@test length(GenerativeAD.Models.latent_score_mean(model,X)) == size(X,2)
	# AAE with VAMP
	k = 4
	pseudoinput_mean = mean(X, dims=ndims(X))
	parameters = (zdim=8, hdim=32, lr=0.001, batchsize=8, activation="swish", nlayers=3, idim=M,
		prior="vamp", pseudoinput_mean=pseudoinput_mean,k=k)
	model = GenerativeAD.Models.aae_constructor(;parameters...)
	model = basic_convergence_test(model)
	@test size(GenerativeAD.Models.reconstruct(model,X)) == size(X)
	@test length(GenerativeAD.Models.reconstruction_score_mean(model,X)) == size(X,2)

	# iris test
	# vanilla AAE
	data = GenerativeAD.load_data("iris")
	include(joinpath(pathof(GenerativeAD), "../../scripts/experiments_tabular/aae.jl"))
	parameters = (zdim=2, hdim=32, lr=0.001, batchsize=32, activation="relu", nlayers=3)
	edited_parameters = GenerativeAD.edit_params(data, parameters)
	training_info, sfs_params = fit(data, edited_parameters)
	save_entries = merge(training_info, (model=nothing,))
	results = map(res->GenerativeAD.experiment(res..., data, ""; save_result=false, save_entries...), sfs_params)
	@info "AAE iris test"
	for result in results
		println(result[:parameters])
		println("validation auc = ", au_roccurve(result[:val_labels], result[:val_scores]))
		println("test auc = ", au_roccurve(result[:tst_labels], result[:tst_scores]))
	end
end


text:
"Square the sum of the first `n` positive integers"
function square_of_sum(n::Integer)::Integer
    # https://www.wolframalpha.com/input/?i=%281+%2B+2+%2B+...+%2B+n%29%5E2
    # (1 + 2 + ... + n)^2 =
    1//4 * n^2 * (n + 1)^2
end

"Sum the squares of the first `n` positive integers"
function sum_of_squares(n::Integer)::Integer
    # https://www.wolframalpha.com/input/?i=1%5E2+%2B+2%5E2+%2B+3%5E2+%2B+...+%2B+n%5E2
    # 1^2 + 2^2 + 3^2 + ... + n^2 =
    1//6 * n * (n + 1) * (2n + 1)
end

"Subtract the sum of squares from square of the sum of the first `n` positive ints"
function difference(n::Integer)::Integer
    # The difference can be further simplified:
    # https://stackoverflow.com/questions/15593428/function-to-calculate-the-difference-between-sum-of-squares-and-square-of-sums/15593631
    -(3n^2 + 2n) * (1 - n^2) // 12
end

max_stars_repo_name: twolodzko/Learning

max_issues_repo_name: twolodzko/Learning

max_forks_repo_name: twolodzko/Learning


text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder

version = v"7.3.0"

# Collection of sources required to build gcc
sources = Any[
    "https://mirrors.kernel.org/gnu/gcc/gcc-7.3.0/gcc-7.3.0.tar.xz" =>
    "832ca6ae04636adbb430e865a1451adf6979ab44ca1c8374f61fba65645ce15c",
]

common_setup = raw"""
    cd $WORKSPACE/srcdir
    cd gcc-7.3.0/
    contrib/download_prerequisites
    export LD_LIBRARY_PATH=/usr/local/lib64:/usr/local/lib:/lib:/usr/glibc-compat/lib:/usr/local/lib:/usr/lib:/opt/x86_64-linux-gnu/lib64:/opt/x86_64-linux-gnu/lib:/opt/x86_64-linux-gnu/x86_64-linux-gnu/lib64:/opt/x86_64-linux-gnu/x86_64-linux-gnu/lib
    cd ..
    mkdir gcc_build
    cd gcc_build
"""

bootstrap_configure = raw"""
../gcc-7.3.0/configure \
    --prefix=/opt/x86_64-linux-gnu \
    --target=x86_64-linux-gnu \
    --host=x86_64-alpine-linux-musl \
    --build=x86_64-alpine-linux-musl \
    --disable-multilib \
    --disable-werror \
    --disable-shared \
    --disable-threads \
    --disable-libatomic \
    --disable-decimal-float \
    --disable-libffi \
    --disable-libgomp \
    --disable-libitm \
    --disable-libmpx \
    --disable-libquadmath \
    --disable-libssp \
    --disable-libsanitizer \
    --without-headers \
    --with-newlib \
    --disable-bootstrap \
    --enable-languages=c \
    --with-sysroot="/opt/x86_64-linux-gnu/x86_64-linux-gnu/sys-root" \
    ${GCC_CONF_ARGS}
"""

main_configure = raw"""
../gcc-7.3.0/configure --prefix=/opt/x86_64-linux-gnu \
    --target=x86_64-linux-gnu \
    --host=x86_64-alpine-linux-musl \
    --build=x86_64-alpine-linux-musl \
    --disable-multilib \
    --disable-werror \
    --enable-languages=c,c++,fortran \
    --with-sysroot=/opt/x86_64-linux-gnu/x86_64-linux-gnu/sys-root \
    --enable-host-shared \
    --enable-threads=posix \
    LD=/opt/super_binutils/bin/
"""

# Bash recipe for building across all platforms
common_suffix = raw"""
make -j${nprocs}
DESTDIR=$WORKSPACE/destdir make install -j20
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = [
    Linux(:x86_64, :musl)
]

# The products that we will ensure are always built
products(prefix) = [
    ExecutableProduct(prefix, "", :xgcc)
]

# Dependencies that must be installed before this package can be built
dependencies = [
#    "https://github.com/staticfloat/GlibcBuilder/releases/download/v2.27-0/build_Glibc.v2.12.2.jl",
    "https://github.com/bicycle1885/ZlibBuilder/releases/download/v1.0.1/build_Zlib.v1.2.11.jl"
]

# Build a boostrap configure
bootstrap_script = """
$common_setup
$bootstrap_configure
$common_suffix
"""
bootstrap_path = joinpath(@__DIR__, "products", "gcc_boostrap.v7.3.0.x86_64-linux-gnu.tar.gz")
if !isfile(bootstrap_path)
    product_hashes = build_tarballs(ARGS, "gcc_boostrap", version, sources, bootstrap_script, platforms, products, dependencies)
    _, bootstrap_hash = product_hashes["x86_64-linux-gnu"]
else
    using SHA: sha256
    bootstrap_hash = open(bootstrap_path) do f
        bytes2hex(sha256(f))
    end
end

push!(sources, Ref("../GlibcBuilder/products/Glibc.v2.12.2.x86_64-linux-gnu.tar.gz" => "7d81df29ec9dfe858244b12a03eb4cc403da84637f458f6e12b6e482dace8d30"))

# Build the tarballs, and possibly a `build.jl` as well.
main_script = """
mkdir -p /opt/x86_64-linux-gnu
tar -C /opt/x86_64-linux-gnu -xvf Glibc.v2.12.2.x86_64-linux-gnu.tar.gz
$common_setup
$main_configure
$common_suffix
"""
build_tarballs(ARGS, "gcc", version, sources, main_script, platforms, products, dependencies)

max_stars_repo_name: Keno/GccBuilder

max_issues_repo_name: Keno/GccBuilder

max_forks_repo_name: Keno/GccBuilder


text:
# isdefined(Base, :__precompile__) && __precompile__()

module VariationalInequality

# package code goes here

using JuMP
import Ipopt

include("model.jl")
include("algorithms.jl")

export  VIPModel,
        VIPData,
        solveVIP,
        correspond,
        clearValues,
        saveSolution,
        gap_function

end # module


text:
"""
`module  JAC.BasicsAZ`
    ... a submodel of JAC that contains methods from the Basics module but with reference to many other modules.
"""
module BasicsAZ

    using Printf,  LinearAlgebra, GaussQuadrature, JAC, ..AlphaVariation, ..AngularMomentum, ..Atomic, ..AtomicState, ..AutoIonization,
                  ..Basics, ..Bsplines, ..Cascade, ..Continuum, ..DecayYield, ..Defaults, ..Dielectronic,
                  ..Einstein, ..FormFactor, ..Hfs, ..HydrogenicIon,
                  ..InteractionStrength, ..IsotopeShift, ..LandeZeeman, ..LSjj, ..ManyElectron, ..MultipolePolarizibility, ..Nuclear,
                  ..PhotoEmission, ..PhotoExcitation, ..PhotoIonization, ..PlasmaShift,
                  ..Radial, ..RadialIntegrals, ..TableStrings

    include("module-BascisAZ-inc-AG.jl")
    include("module-BascisAZ-inc-compute.jl")
    include("module-BascisAZ-inc-generate.jl")
    include("module-BascisAZ-inc-HP.jl")
    include("module-BascisAZ-inc-perform.jl")
    include("module-BascisAZ-inc-QZ.jl")

end # module


text:
module SpectralSolver

using LinearAlgebra

function Fₛ(s,l,m)
    if l < min(abs(s),abs(m))
        return 0.0
    else
        num1 = (l+1)^2 - m^2
        num2 = (l+1)^2 - s^2
        denom = (2*l+3)*(2*l+1)*(l+1)^2
        return sqrt(num1*num2/denom)
    end
end

function Gₛ(s,l,m)
    if l == zero(l)
        return 0.
    else
        num1 = l^2 - m^2
        num2 = l^2 - s^2
        denom = (4*l^2 - 1)*l^2
        return sqrt(num1*num2/denom)
    end
end

function Hₛ(s,l,m)
    if (l == 0) & (s == 0)
        return 0.
    else
        return (-m*s)/(l*(l+1))
    end
end

Aₛ(s,l,m) = Fₛ(s,l,m)*Fₛ(s,l+1,m)
Dₛ(s,l,m) = Fₛ(s,l,m)*(Hₛ(s,l+1,m) + Hₛ(s,l,m))
Bₛ(s,l,m) = Fₛ(s,l,m)*Gₛ(s,l+1,m) + Gₛ(s,l,m)*Fₛ(s,l-1,m) + Hₛ(s,l,m)^2
Eₛ(s,l,m) = Gₛ(s,l,m)*(Hₛ(s,l-1,m) + Hₛ(s,l,m))
Cₛ(s,l,m) = Gₛ(s,l,m)*Gₛ(s,l-1,m)

function M(l,lʼ,s,m,c)
    if lʼ == l-2
        return (-c^2)*(Aₛ(s,lʼ,m))
    elseif lʼ == l-1
        return (-c^2)*(Dₛ(s,lʼ,m)) + 2*c*s*Fₛ(s,lʼ,m)
    elseif lʼ == l
        return lʼ*(lʼ + 1) - s*(s+1) - (c^2)*Bₛ(s,lʼ,m) + 2*c*s*Hₛ(s,lʼ,m)
    elseif lʼ == l+1
        return (-c^2)*(Eₛ(s,lʼ,m)) + 2*c*s*Gₛ(s,lʼ,m)
    elseif lʼ == l+2
        return (-c^2)*(Cₛ(s,lʼ,m))
    else
        return 0.0*im
    end
end

lmin(s,m) = max(abs(m),abs(s))

function ConstructM(lmax,s,m,c)
    lmins = lmin(m,s)
    N = lmax - lmin(m,s) + 1
    Mₑ = zeros(Complex{Float64},N,N)
    for i = 1:N, j = 1:N
        l = i+lmins-1; lʼ = j+lmins-1;
        if abs(l-lʼ) <= 2
            Mₑ[i,j] = M(l,lʼ,s,m,c)
        end
    end
    Mₑ
end

Mₑ = ConstructM(10,2,-1,0.01+0.01*im)
vals, vecs = eigen(Mₑ)

function ComputeAₗₘ(s,m,c, Aₗₘprev,lmax)
    Mₑ = ConstructM(lmax,s,m,c)
    vals, vecs = eigen(Mₑ)
    error, j = findmin(map(Anew -> abs2(Anew - Aₗₘprev), vals))
    vals[j], vecs[:,j], error
end

export ComputeAₗₘ
end

max_stars_repo_name: Potatoasad/AttemptAtQNM

max_issues_repo_name: Potatoasad/AttemptAtQNM

max_forks_repo_name: Potatoasad/AttemptAtQNM


text:
using PkgBenchmark

mkconfig(; kwargs...) =
    BenchmarkConfig(
        env = Dict(
            "JULIA_NUM_THREADS" => "1",
            "OMP_NUM_THREADS" => "1",
        );
        kwargs...
    )

progressoptions =
    if lowercase(get(ENV, "CI", "false")) == "true"
        (dt = 60 * 9.0,)
    else
        NamedTuple()
    end

group_target = benchmarkpkg(
    dirname(@__DIR__),
    mkconfig(),
    progressoptions = progressoptions,
    resultfile = joinpath(@__DIR__, "result-target.json"),
)

group_baseline = benchmarkpkg(
    dirname(@__DIR__),
    mkconfig(id = "master"),
    progressoptions = progressoptions,
    resultfile = joinpath(@__DIR__, "result-baseline.json"),
)

judgement = judge(group_target, group_baseline)


text:
using BinaryBuilder
using GitHub
using Pkg
using Random
using ghr_jll

const BB = BinaryBuilder
const GH = GitHub

function extract_targets!(args)
    targets = filter(arg -> !startswith(arg, "--"), args)
    # Remove collected targets from ARGS as we handle them directly
    filter!(arg -> startswith(arg, "--"), args)
    return targets
end

function save(output, name = "build_output.jl")
    # Indent levels
    s¹ = "    "
    s² = s¹ * s¹
    s³ = s² * s¹
    s⁴ = s² * s²
    s⁵ = s³ * s²

    open(name, "w") do io
        write(io, "using BinaryBuilder\n")
        write(io, "using Pkg\n\n")
        write(io, "using Base: SHA1\n\n")
        write(io, "output = Dict(\n")
        for version in keys(output)
            write(io, s¹, repr(version), " => Dict(\n")
            platforms_products = output[version]
            for platform in keys(platforms_products)
                p = platforms_products[platform]
                write(io, s², repr(platform), " => (\n")
                write(io, s³, repr(p[1]), ",\n")
                write(io, s³, repr(p[2]), ",\n")
                write(io, s³, repr(p[3]), ",\n")
                write(io, s³, "Dict(\n")
                products = p[4]
                for product in keys(products)
                    # The `LibraryProduct` constructor does not match its string
                    # representation
                    if product isa BinaryBuilder.LibraryProduct
                        write(io, s⁴, "LibraryProduct(\n")
                        write(io, s⁵, repr(product.libnames), ", ")
                        write(io, repr(product.variable_name), ", ")
                        write(io, repr(product.dir_paths), ";\n")
                        write(io, s⁵, "dont_dlopen = ", repr(product.dont_dlopen), ", ")
                        write(io, "dlopen_flags = ", repr(product.dlopen_flags), "\n")
                        write(io, s⁴, ")")
                    else
                        write(io, s⁴, repr(product))
                    end
                    write(io, " => Dict(\n")
                    for pair in products[product]
                        write(io, s⁵, repr(pair), ",\n")
                    end
                    write(io, s⁴, "),\n")
                end
                write(io, s³, "),\n")
                write(io, s², "),\n")
            end
            write(io, s¹, "),\n")
        end
        write(io, ")")
    end
end

function build_jlls(ARGS, name, output, dependencies)
    args = copy(ARGS)
    verbose = BB.check_flag!(args, "--verbose")
    gh_account = last(BB.extract_flag!(args, "--gh-account", "JuliaBinaryWrappers"))

    jll_repo = "$(gh_account)/$(name)_jll.jl"
    jll_path = joinpath(Pkg.devdir(), "$(name)_jll")

    BB.init_jll_package(name, jll_path, jll_repo)

    versions = (sort! ∘ collect ∘ keys)(output)
    for version in versions
        platforms_products = output[version]
        build_version = get_next_version(jll_repo, version)
        build_tag = "v$(build_version)"
        url = "https://github.com/$(jll_repo)/releases/download/$(build_tag)"

        if verbose
            message = (
                "Committing and pushing $(name)_jll.jl wrapper code version " *
                "$(build_version)..."
            )
            @info(message)
        end
        BB.build_jll_package(
            name, build_version, jll_path, platforms_products, dependencies, url;
            verbose = verbose
        )
        BB.push_jll_package(
            name, build_version; code_dir = jll_path, deploy_repo = jll_repo
        )

        # Copy current version tarballs to a temporary path to prevent tarball mixing
        name_version = "$(name).v$(version)"
        products_path = joinpath(pwd(), "products")
        tmp_path = joinpath(products_path, randstring())
        mkpath(tmp_path)
        for tarball in readdir(products_path)
            if startswith(tarball, name_version)
                src = joinpath(products_path, tarball)
                dst = joinpath(tmp_path, tarball)
                cp(src, dst; force = true)
            end
        end

        if verbose
            @info("Deploying tarballs to release $(build_tag) on $(jll_repo) via `ghr`...")
        end
        BB.upload_to_github_releases(jll_repo, build_tag, tmp_path; verbose = verbose)

        # Clean up
        rm(tmp_path; recursive = true, force = true)
    end
end

function get_next_version(repo, version)
    # If version already has a build_number, try using it
    build_number = UInt(0)
    if version.build != ()
        build_number = first(version.build)
    end

    # Download tag list from GitHub
    refs = first(GH.gh_get_paged_json(GH.DEFAULT_API, "/repos/$(repo)/git/refs"))
    filter!(r -> occursin(r"tags/v", r["ref"]), refs)

    # Collect the versions that match ours
    versions = VersionNumber[]
    for r in refs
        v = VersionNumber(last(rsplit(r["ref"], '/'; limit = 2)))
        version_matches = (
            v.major == version.major && v.minor == version.minor &&
            v.patch == version.patch && v.build isa Tuple{UInt}
        )
        if version_matches
            push!(versions, v)
        end
    end

    # Our build number must be larger that the maximum already deployed
    if !isempty(versions)
        last_build_number = first(maximum(versions).build)
        if build_number ≤ last_build_number
            build_number = last_build_number + UInt(1)
        end
    end

    return VersionNumber(
        version.major, version.minor, version.patch, version.prerelease, (build_number,)
    )
end

max_stars_repo_name: pabloferz/LAMMPSBuilder

max_issues_repo_name: pabloferz/LAMMPSBuilder

max_forks_repo_name: pabloferz/LAMMPSBuilder


text:
"""
    smallest_subset4sum(array::Array{Int64,1}, target::Int64)

Returns length of smallest subarray-unit with sum greater than the target.

# Arguments
- `array::Array{Int64,1}`: Unsorted array
- `target::Int64`: Target number to reach

# Examples
```julia-repl
julia> import ClassicAlgorithmsCollections
julia> arr =  [10, 12, 11, 9, 13, 14, 17, 18, 15]
julia> ClassicAlgorithmsCollections.smallest_subarray4sum(arr)
4
```
"""
function smallest_subset4sum(array::Array{Int64,1}, target::Int64)
    n = length(array)
    min_length = n

    for i in 1:n
        current_sum = array[i]
        if current_sum > target
            return 1
        end
        for j in (i+1):n
            current_sum += array[j]
            if current_sum > target && (j - i + 1) < min_length
                min_length = (j - i + 1)
            end
        end
    end
    if min_length != n
        return min_length
    else
        return nothing
    end
end

max_stars_repo_name: Anselmoo/ClassicAlgorithmsCollections

max_issues_repo_name: Anselmoo/ClassicAlgorithmsCollections

max_forks_repo_name: Anselmoo/ClassicAlgorithmsCollections


text:
# Think about updating the reference trial writing
# to write reference trials to directories with the name of the git commit
# that the trial corresponds to.
# I.e. if the trial was run under commit f85e8e7, then within the top-level
# "ref" directory, there will be a directories named after commits
# Will also need to think about best ways to parse these directories
# in some sort of chronological order, or having a bunch of commit id labeled
# directories could get a bit unwieldy

using JLD2

"""
```
write_ref_trial(trial, trial_name)
```

Write a reference trial to a JLD file, to act as the standard that new trials are
benchmarked against.

### Arguments
- `trial::BenchmarkTools.Trial`: The trial object that is being written.
- `trial_name::String`: The name of the trial being written.

"""
function write_ref_trial(trial::BenchmarkTools.Trial, trial_name::String)
    filepath = "../reference/$trial_name.jld"
    if isfile(filepath)
        println("There is already a reference file in $filepath. Do you want to overwrite it? (y/n)")
        response = readline(STDIN)
        if response == "y"
            rm(filepath)
            op = "Overwrote"
        elseif response == "n"
            println("File will be left as is. Aborting write.")
            return nothing
        else
            throw("Invalid response, must be y or n")
        end
    else
        op = "Wrote"
    end

    d = Dict{Symbol, Union{Vector, Int}}()
    d[:times] = trial.times
    d[:gctimes] = trial.gctimes
    d[:memory] = trial.memory
    d[:allocs] = trial.allocs

    save(filepath, trial_name, d)

    println("$op "*filepath)
end

function write_ref_trial_group(group::Dict{Symbol, Dict},
                               trial_names::Vector{Symbol},
                               group_name::String)
    filepath = "../reference/$group_name.jld"
    if isfile(filepath)
        println("There is already a reference file in $filepath. Do you want to overwrite it? (y/n)")
        response = readline(STDIN)
        if response == "y"
            rm(filepath)
            op = "Overwrote"
        elseif response == "n"
            println("File will be left as is. Aborting write.")
            return nothing
        else
            throw("Invalid response, must be y or n")
        end
    else
        op = "Wrote"
    end

    save(filepath, group_name, group)
    println("$op "*filepath)
end


text:
"""
zero kernel function
"""
function zero_kernel(
    hyperparameters::Vector{<:Real},
    δ::Real;
    dorder::Vector{<:Integer}=zeros(Int64, length(hyperparameters) + 2))

    return 0

end

return zero_kernel, 0  # the function handle and the number of kernel hyperparameters

max_stars_repo_name: christiangil/GPLinearODEMaker

max_issues_repo_name: christiangil/GPLinearODEMaker

max_forks_repo_name: christiangil/GPLinearODEMaker


text:
"""
Use jld macros @save and @load... Works unreliably.
https://stackoverflow.com/questions/40970732/julia-save-workspace

User points out that these macros are not meant for within function scope.
"""

using QuNet
using LightGraphs, SimpleWeightedGraphs, GraphPlot, MetaGraphs, Plots, Colors, Statistics
using LaTeXStrings
# using DelimitedFiles
using JLD

"""
Plot the performance statistics of greedy_multi_path for 1 random end-user pair
over a range of graph sizes and for different numbers of max_path (1,2,3) (ie. the
maximum number of paths that can be purified by a given end-user pair.)
"""

generate_new_data = false
if generate_new_data == true

    # Params
    num_trials = 1000::Int64
    min_size = 2::Int64
    max_size = 20::Int64

    @assert min_size < max_size

    perf_data1 = []
    perf_data2 = []
    perf_data3 = []
    perf_data4 = []

    size_list = collect(min_size:1:max_size)
    for (j, data_array) in enumerate([perf_data1, perf_data2, perf_data3, perf_data4])
        println("Collecting for max_paths: $j")
        for i in size_list
            println("Collecting for gridsize: $i")
            # Generate ixi graph:
            net = GridNetwork(i, i)

            # Collect performance statistics
            performance, dummy, dummy, dummy = net_performance(net, num_trials, 1, max_paths=j)
            push!(data_array, performance)
        end
    end

    # Get values for x axis
    x = collect(min_size:1:max_size)

    # Save variables to julia data file
    @save "test/maxpaths.jdl"
end

# Load variables from data file
@load "test/maxpaths.jdl"

# debug
println("Environment check")
println(num_trials)

# Extract data from performance data
loss_arr1 = collect(map(x->x["loss"], perf_data1))
z_arr1 = collect(map(x->x["Z"], perf_data1))
loss_arr2 = collect(map(x->x["loss"], perf_data2))
z_arr2 = collect(map(x->x["Z"], perf_data2))
loss_arr3 = collect(map(x->x["loss"], perf_data3))
z_arr3 = collect(map(x->x["Z"], perf_data3))
loss_arr4 = collect(map(x->x["loss"], perf_data4))
z_arr4 = collect(map(x->x["Z"], perf_data4))

# Plot
plot(x, loss_arr1, linewidth=2, label=L"$\eta_1$", xlims=(0, max_size), color = :red, legend=:right)
plot!(x, z_arr1, linewidth=2, label=L"$F_1$", linestyle=:dash, color =:red)
plot!(x, loss_arr2, linewidth=2, label=L"$\eta_2$", color =:blue)
plot!(x, z_arr2, linewidth=2, label=L"$F_2$", linestyle=:dash, color =:blue)
plot!(x, loss_arr3, linewidth=2, label=L"$\eta_3$", color =:green)
plot!(x, z_arr3, linewidth=2, label=L"$F_3$", linestyle=:dash, color =:green)
plot!(x, loss_arr4, linewidth=2, label=L"$\eta_3$", color =:purple)
plot!(x, z_arr4, linewidth=2, label=L"$F_3$", linestyle=:dash, color =:purple)
xaxis!(L"$\textrm{Grid Size}$")

# make analytic function
function M(n)
    return 2/3 * n * (n^2 - 1) / (2 * n^2 - 1)
end

function E(n)
    return 10^(-M(n) / 10)
end

function F(n)
    return (10^(-M(n) / 10) + 1)/2
end

x = collect(min_size:1:max_size)
ave_e = E.(x)
ave_f = F.(x)

plot!(x, ave_e, linewidth=2, label=L"$\textrm{ave } \eta$", color =:orange)
plot!(x, ave_f, linewidth=2, label=L"$\textrm{ave } F$", linestyle=:dash, color =:orange)

savefig("cost_maxpaths.png")
savefig("cost_maxpaths.pdf")

# Plot
plot(x[], loss_arr1, linewidth=2, label=L"$\eta_1$", xlims=(0, max_size), yaxis=:log, color = :red, legend=:left)
plot!(x, loss_arr2, linewidth=2, label=L"$\eta_2$", yaxis=:log, color =:blue)
plot!(x[3:end], loss_arr3[3:end], linewidth=2, label=L"$\eta_3$", yaxis=:log, color =:green)
plot!(x[3:end], loss_arr4[3:end], linewidth=2, label=L"$\eta_4$", yaxis=:log, color =:purple)

z_arr1 = z_arr1 .- 0.5
z_arr2 = z_arr2 .- 0.5
z_arr3 = z_arr3 .- 0.5
z_arr4 = z_arr4 .- 0.5

plot!(x[3:end], z_arr1[3:end], linewidth=2, label=L"$F_1 - 0.5$", linestyle=:dash, color =:red, yaxis=:log)
plot!(x[3:end], z_arr2[3:end], linewidth=2, label=L"$F_2 - 0.5$", linestyle=:dash, color =:blue, yaxis=:log)
plot!(x[3:end], z_arr3[3:end], linewidth=2, label=L"$F_3 - 0.5$", linestyle=:dash, color =:green, yaxis=:log)
plot!(x[3:end], z_arr4[3:end], linewidth=2, label=L"$F_4 - 0.5$", linestyle=:dash, color =:purple, yaxis=:log)
xaxis!(L"$\textrm{Grid Size}$")

savefig("cost_maxpathslog.png")
savefig("cost_maxpathslog.pdf")

max_stars_repo_name: FalafelGood/QuNet

max_issues_repo_name: FalafelGood/QuNet

max_forks_repo_name: FalafelGood/QuNet


text:
using PkgSRA, Plots, LinearAlgebra, Statistics
using BSON: @load
pyplot()

#####
##### Load variables
#####
include("paper_settings.jl")
this_dat_name = DAT_FOLDERNAME*"TMP_dat_flowchart_inkscape_"

# Uncontrolled data
fname = this_dat_name*"uncontrolled_lorenz.bson"
@load fname dat_raw numerical_grad_raw

# Original ODE and controller variables
fname = this_dat_name*"ode_vars.bson";
@load fname dat dat_grad true_grad U_true

# Bayesian variables 1: before control
fname = this_dat_name*"naive_vars_sindy.bson";
@load fname sindy_dat_unctr sindy_unctr

# Bayesian variables 1: before control
fname = this_dat_name*"naive_vars_bayes.bson";
@load fname dat_unctr sample_trajectory_noise residual accepted_ind ctr_guess

# SINDY variables 2: after control
fname = this_dat_name*"ctr_vars_sindy.bson";
@load fname sindy_sub sindy_dat_ctr sindy_grad_ctr

# Bayesian variables 2: after control
fname = this_dat_name*"ctr_vars_turing.bson";
@load fname dat_ctr sample_trajectory_noise_ctr sample_gradients_ctr turing_dat_ctr ctr_guess2

#####
##### Plot everything
#####
this_fig_name = FIGURE_FOLDERNAME*"TMP_fig_flowchart_inkscape_"
plot_opt = Dict(:titlefontsize=>48,
        :xticks=>false, :yticks=>false, :zticks=>false,
        :legend=>false, :fontfamily=>:serif)

## 0: Raw, uncontrolled system
plot_raw = plot3d(dat_raw[1, :], dat_raw[2, :], dat_raw[3, :],
        color=COLOR_DICT["data_uncontrolled"], lw=4; plot_opt...)
        title!("Unforced Lorenz")
fname = this_fig_name * "uncontrolled.png";
savefig(plot_raw, fname)


text:
error_green_formula(SL,DL,γ₀u,γ₁u,u,σ)                      = σ*u + SL*γ₁u  - DL*γ₀u
error_derivative_green_formula(ADL,H,γ₀u,γ₁u,un,σ)          = σ*un + ADL*γ₁u - H*γ₀u
error_interior_green_identity(SL,DL,γ₀u,γ₁u)                = error_green_formula(SL,DL,γ₀u,γ₁u,γ₀u,-1/2)
error_interior_derivative_green_identity(ADL,H,γ₀u,γ₁u)     = error_derivative_green_formula(ADL,H,γ₀u,γ₁u,γ₁u,-1/2)
error_exterior_green_identity(SL,DL,γ₀u,γ₁u)                = error_green_formula(SL,DL,γ₀u,γ₁u,γ₀u,1/2)
error_exterior_derivative_green_identity(ADL,H,γ₀u,γ₁u)     = error_derivative_green_formula(ADL,H,γ₀u,γ₁u,γ₁u,1/2)

max_stars_repo_name: WaveProp/WaveProp

max_issues_repo_name: WaveProp/WaveProp

max_forks_repo_name: WaveProp/WaveProp


text:
function en(an,anm1,Ln)
	return sqrt((1-abs2(an))/((one(an)-abs2(anm1))*(one(an)-abs2(Ln))))
end

function ln(z::Vector{Complex128},an::Complex128,anm1::Complex128,phinm1::Vector{Complex128},phinm1star::Vector{Complex128})
    f1 = [((z[i]-anm1)/(one(Complex128)-conj(an)*z[i]))*conj(phinm1[i]) for i=1:length(z)]
    f2 = [((one(Complex128)-conj(anm1)*z[i])/(one(Complex128)-conj(an)*z[i]))*conj(phinm1star[i]) for i=1:length(z)]
    num = sum(f1)
    den = sum(f2)
    return -num/den
end

function phiNext(z::Vector{Complex128},an::Complex128,anm1::Complex128,phinm1::Vector{Complex128},phinm1star::Vector{Complex128})
    phin = zeros(Complex128,length(phinm1))
    phinstar = zeros(Complex128,length(phinm1))
    L = ln(z,an,anm1,phinm1,phinm1star)
    c = en(an,anm1,L)
    for k=1:length(z)
        q = (z[k]-anm1)
        p = (one(Complex128)-conj(anm1)*z[k])
        C = c/(z[k]-an)
        PHI = p*phinm1[k]
        PHISTAR = q*phinm1star[k]
        phin[k] = C*(PHI+conj(L)*PHISTAR)
        phinstar[k] = C*(L*PHI+PHISTAR)
    end
    return (phin, phinstar)
end

function genZ(N::Int)
	return exp((im*2*pi/N)*[0:N-1])
end

function genBasis(N::Int,a::Vector{Complex128})
    z = genZ(N)
    M = length(a)+1
    a = cat(1,zero(Complex128),a)
    B = ones(Complex128,(N,M))/sqrt(N)
    BstarOld = ones(Complex128,N)/sqrt(N)
    for k=2:M
        (B[:,k],BstarNew) = phiNext(z,a[k],a[k-1],B[:,k-1],BstarOld)
        BstarOld = BstarNew
    end
    return B[:,2:end]
end

function genWBasis(N::Int,a::Vector{Complex128})
    z = genZ(N)
    M = length(a)+1
    a = cat(1,zero(Complex128),a)
    B = ones(Complex128,(N,M))/sqrt(N)
    BstarOld = ones(Complex128,N)/sqrt(N)
    for k=2:M
        (B[:,k],BstarNew) = phiNext(z,a[k],a[k-1],B[:,k-1],BstarOld)
        BstarOld = BstarNew
    end
    return fftshift(B[:,2:end],1)
end

function genTBasis(N::Int,a::Vector{Complex128})
	B = genBasis(N,a)
	T = fft(B,1)/sqrt(N)
    T = flipud(T)
    return T
end

function genA(N::Int)
	return (rand(N)).*exp(im*2*pi*rand(N))
end

function compress2D(ay,ax,D)
	(y,x) = size(D)
	W = fftshift(genBasis(y,ay),1)
	T = genTBasis(x,ax)
	return ((W'*(D*T)),W,T)
end
fftshift
function uncompress2D(X,W,T)
	return W*X*T'
end

function lsq2DTMobj(beta::Vector,grad::Vector,Nx::Int,Ny::Int,D::Array{Complex128,2})
    beta1 = beta[1:(2*Ny)]
    beta1c = [beta1[i]*exp(im*2*pi*beta1[i+1]) for i=1:2:length(beta1)]
    beta2 = beta[(2*Ny+1):end]
    beta2c = [beta2[i]*exp(im*2*pi*beta2[i+1]) for i=1:2:length(beta2)]
    (X,W,T) = compress2D(beta1c,beta2c,D)
    Dhat = uncompress2D(X,W,T)
    return norm(D-Dhat)
end

function lsq2DTMgenBounds(Nx::Int,Ny::Int,offset)
    betaLB = zeros(Float64,2*Nx+2*Ny)
    betaUB = ones(Float64,2*Nx+2*Ny)
    betaUB[1:2:end] -= offset
    return (betaLB,betaUB)
end

function lsq2DTMgenStart(Nx::Int,Ny::Int,offset::Float64)
    Nt = 2*Nx+2*Ny
    Nr = Nx+Ny
    betaGuess = zeros(Float64,Nt)
    betaGuess[1:2:end] = (one(Float64)-offset)*rand(Nr)
    betaGuess[2:2:end] = rand(Nr)
    return betaGuess
end

function lsq2DTMopt(Nx::Int,Ny::Int,D::Array{Complex128,2},gmethod,lmethod,ltol,gWallClockTime)
    #This function assumes that NLopt is in the namespace
    Nt = 2*Nx+2*Ny
    offset = 0.001
    opt = Opt(gmethod,Nt)
    lopt = Opt(lmethod,Nt)
    ftol_rel!(lopt, ltol) #local stopping criteria
    maxtime!(opt,gWallClockTime) #global stopping critera
    (betalb,betaub) = lsq2DTMgenBounds(Nx,Ny,offset)
    beta0 = lsq2DTMgenStart(Nx,Ny,offset)
    local_optimizer!(opt,lopt)
    min_objective!(opt,(x,g)->lsq2DTMobj(x,g,Nx,Ny,D))
    lower_bounds!(opt,betalb)
    upper_bounds!(opt,betaub)
    (rss,betaopt,exitcode) = optimize(opt,beta0)
    betaopty = betaopt[1:(2*Ny)]
    betaoptx = betaopt[(2*Ny+1):end]
    betay = [betaopty[i]*exp(im*2*pi*betaopty[i+1]) for i=1:2:length(betaopty)]
    betax = [betaoptx[i]*exp(im*2*pi*betaoptx[i+1]) for i=1:2:length(betaoptx)]
    Wopt = genBasis(size(D,1),betay)
    Topt = genTBasis(size(D,2),betax)
    #Xopt = (W'*(D*T))
    return (betay,betax,Wopt,Topt,rss,exitcode)
end

function phiNextNoStar(z::Vector{Complex128},an::Complex128,anm1::Complex128,phinm1::Vector{Complex128},phinm1star::Vector{Complex128})
    phin = zeros(Complex128,length(phinm1))
    L = ln(z,an,anm1,phinm1,phinm1star)
    c = en(an,anm1,L)
    for k=1:length(z)
        q = (z[k]-anm1)
        p = (one(Complex128)-conj(anm1)*z[k])
        C = c/(z[k]-an)
        PHI = p*phinm1[k]
        PHISTAR = q*phinm1star[k]
        phin[k] = C*(PHI+conj(L)*PHISTAR)
    end
    return phin
end

function unitDiskGrid(n::Int,delta::Float64)
    N = nextpow2(n)
    M = N>>1
    v = linspace(delta,one(Float64)-delta,M)
    Zjk = zeros(Complex128,N,N)
    for j=1:(M)
            #bottom, running counter-clockwise
            Zjk[(M+j),((M-(j-1)):(M+1+(j-1)))] = v[j]*exp(im*linspace(5*pi/4,7*pi/4,2*j))
            #top, running clockwise
            Zjk[(M-j+1),((M-(j-1)):(M+1+(j-1)))] = v[j]*exp(im*linspace(3*pi/4,1*pi/4,2*j))
    end
    for j=1:(M)
            #right, running clockwise
            Zjk[((M-(j-1)):(M+1+(j-1))),M+j] = v[j]*exp(im*linspace(pi/4,-pi/4,2*j)) #right
            #left, running counter-clockwise
            Zjk[((M-(j-1)):(M+1+(j-1))),M-j+1] = v[j]*exp(im*linspace(3*pi/4,5*pi/4,2*j)) #left
    end
    return (reshape(Zjk,N^2))
end

function genOMPBasis(N::Int,Z::Vector{Complex128},anm1::Complex128,vB::Vector,vBstar::Vector)
    z = genZ(N)
    M = length(Z)
    B = zeros(Complex128,(N,M))
    for k=1:M
        B[:,k] = phiNextNoStar(z,Z[k],anm1,vB,vBstar)
    end
    return B
end

function genOMPTBasis(N::Int,Z::Vector{Complex128},anm1::Complex128,vB::Vector,vBstar::Vector)
    B = genOMPBasis(N,Z,anm1,vB,vBstar)
    T = fft(B,1)/sqrt(N)
    T = flipud(T) #flip to be causal, instead of acausal
    return T
end

function indmaxmat(M)
    (x,y) = size(M)
    linind = indmax(M)
    j = linind%x
    k = fld(linind,y)+1
    return (j,k)
end

function indminmat(M)
    (x,y) = size(M)
    linind = indmin(M)
    j = linind%x
    k = fld(linind,y)+1
    return (j,k)
end

function lsq2DTMOMPobjc(Z::Vector,D::Array{Complex128,2},pym1::Complex128,pxm1::Complex128,vBY::Vector,vBstarY::Vector,vBX::Vector,vBstarX::Vector)
    N = length(Z)
    by = genOMPBasis(size(D,1),Z,pym1,vBY,vBstarY)
    by = fftshift(by,1)
    bx = genOMPTBasis(size(D,2),Z,pxm1,vBX,vBstarX)
    O = (by'*(D*bx))
    (j,k) = indmaxmat(abs2(O))
    (bxopt,bxstaropt) = phiNext(genZ(size(D,2)),Z[k],pxm1,vBX,vBstarX)
    (byopt,bystaropt) = phiNext(genZ(size(D,1)),Z[j],pym1,vBY,vBstarY)
    return (j,k,byopt,bystaropt,bxopt,bxstaropt,O[j,k],O)
end

function lsq2DTMOMPobjls(Z::Vector,D::Array{Complex128,2},ay::Vector,ax::Vector)
    N = length(Z)
    O = zeros(Float64,(N,N))
    (y,x) = size(D)
    for j=1:N
        aynew = [ay; Z[j]]
        BY = fftshift(genBasis(y,aynew),1)
        for k=1:N
            axnew = [ax; Z[k]]
            BX = genTBasis(x,axnew)
            C = BY'*(D*BX)
            O[j,k] = norm(C)
        end
    end
    (jopt,kopt) = indmaxmat(O)
    return (Z[jopt],Z[kopt])
end

function lsq2DTMOMPresid(ay::Vector,ax::Vector,D::Array{Complex128,2})
    #return the residual of an omp step
    BY = fftshift(genBasis(size(D,1),ay),1)
    BX = genTBasis(size(D,2),ax)
    C = BY'*(D*BX)
    Dhat = BY*(C*BX')
    return D-Dhat
end

function lsq2DTMOMPiterate(Niter::Int,Ngrid::Int,D::Array{Complex128,2})
    polesx = zeros(Complex128,Niter+1)
    polesy = zeros(Complex128,Niter+1)
    (y,x) = size(D)
    Z = unitDiskGrid(nextpow2(int(floor(sqrt(Ngrid)))),0.001)
    for i=2:(Niter+1)
        (polesy[i],polesx[i]) = lsq2DTMOMPobjls(Z,D,polesy[2:i],polesy[2:i])
    end
    return (polesy[2:end],polesx[2:end])
end

function separated2DTMgenBounds(N::Int,offset::Float64)
    betaLB = zeros(Float64,2*N)
    betaUB = ones(Float64,2*N)
    betaUB[1:2:end] -= offset
    return (betaLB,betaUB)
end

function separated2DTMgenStart(N::Int,offset::Float64)
    Nt = 2*N
    betaGuess = zeros(Float64,Nt)
    betaGuess[1:2:end] = (one(Float64)-offset)*rand(N)
    betaGuess[2:2:end] = rand(N)
    return betaGuess
end

function separated2DTMobj(beta::Vector,g::Vector,D::Vector{Complex128},f::Function)
    betac = [beta[i]*exp(im*2*pi*beta[i+1]) for i=1:2:length(beta)]
    B = f(size(D,1),betac)
    return norm(D-B*(B'*D))
end

function separated2DTMopt(N::Int,D::Vector{Complex128},basisGen::Function,gmethod,lmethod,ltol,gWallClockTime)
    #This function assumes that NLopt is in the namespace
    offset = 0.001
    opt = Opt(gmethod,2*N)
    lopt = Opt(lmethod,2*N)
    ftol_rel!(lopt, ltol) #local stopping criteria
    maxtime!(opt,gWallClockTime) #global stopping critera
    (betalb,betaub) = separated2DTMgenBounds(N,offset)
    beta0 = separated2DTMgenStart(N,offset)
    local_optimizer!(opt,lopt)
    min_objective!(opt,(x,g)->separated2DTMobj(x,g,D,basisGen))
    lower_bounds!(opt,betalb)
    upper_bounds!(opt,betaub)
    (rss,betaopt,exitcode) = optimize(opt,beta0)
    betaoptc = [betaopt[i]*exp(im*2*pi*betaopt[i+1]) for i=1:2:length(betaopt)]
    return betaoptc
end

function separated2DTM(Npoles::Int,Nsing::Int,D::Array{Complex128,2},TperOpt)
    (u,s,v) = svd(D)
    vp = v[:,1:Nsing]
    sp = s[1:Nsing]
    up  = u[:,1:Nsing]
    betaoptv = zeros(Complex128,(Nsing,Npoles))
    betaoptu = zeros(Complex128,(Nsing,Npoles))
    for k=1:Nsing
        betaoptv[k,:] = separated2DTMopt(Npoles,vp[:,k],genTBasis,:G_MLSL_LDS,:LN_SBPLX,1e-6,TperOpt)
    end
    for k=1:Nsing
        betaoptu[k,:] = separated2DTMopt(Npoles,up[:,k],genWBasis,:G_MLSL_LDS,:LN_SBPLX,1e-6,TperOpt)
    end
    return (betaoptv,betaoptu)
end


text:
module FormFactor

using SpecialFunctions
using HypergeometricFunctions
using WignerFamilies
using FLoops

include("extra_form.jl")

export mform,mfun

end


text:
module WebAssemblyText

using Core: Compiler, SSAValue, SlotNumber, TypedSlot, GotoNode, CodeInfo, NewvarNode
using Core.Compiler: GotoIfNot, PartialStruct
using Core: ReturnNode

"""
Evalscope is a dummy module where all evaluated expressions live.
"""
module Evalscope end

mutable struct BlockInfo
    gotos::Dict{Int,Int}
    parents::Array{Array{Int,1},1}
end

include("structure.jl")
include("itemtype.jl")
include("utils.jl")
include("translate.jl")
include("controlflow.jl")
include("stringify.jl")
include("builtins.jl")
include("imports.jl")
include("debug.jl")
include("parse.jl")
include("jl2wat.jl")

export jl2wat, jlstring2wat, @code_wat

end


text:
facts("Testing random agent selection") do
    agents = TcrCore.setupAgentsWithFixedAccuracy(1, 100)
    onlyAgent = agents[1]
    @fact TcrCore.selectRandomAgent(agents) --> onlyAgent
end

max_stars_repo_name: mlodykaz/tcr-simulation

max_issues_repo_name: mlodykaz/tcr-simulation

max_forks_repo_name: mlodykaz/tcr-simulation


text:
module TestLibrary
include("preamble.jl")
using Dates
using BangBang: push!!

@testset "Cat" begin
    # Inner transducer is stateful:
    @testset for xs in iterator_variants([1:5, 1:2])
        xf = DropWhile(x -> x < 3)
        @test xs |> Cat() |> xf |> collect == collect(xf, vcat(xs...))
    end

    # Inner transducer is stateful and requires complete:
    @testset for xs in iterator_variants([1:5, 1:2])
        xf = opcompose(PartitionBy(isequal(3)), Map(copy))
        @test xs |> Cat() |> xf |> collect == collect(xf, vcat(xs...))
    end
end

@testset "MapCat" begin
    # https://clojuredocs.org/clojure.core/mapcat
    # (transduce (mapcat reverse) conj [[3 2 1 0] [6 5 4] [9 8 7]])
    @testset for xs in iterator_variants(
            [[3, 2, 1, 0], [6, 5, 4], [9, 8, 7]]
            )
        @test collect(MapCat(reverse), xs) == 0:9
        @test xs |> MapCat(reverse) |> Map(inc) |> collect == 1:10
        @test xs |> MapCat(x -> reverse(x) |> Map(inc)) |> collect == 1:10
    end
end

@testset "PartitionBy" begin
    # https://clojuredocs.org/clojure.core/partition-by#example-542692c7c026201cdc3269da
    xf = opcompose(PartitionBy(isequal(3)), Map(copy))
    @testset for xs in iterator_variants(1:5)
        @test collect(xf, xs) == [[1, 2], [3], [4, 5]]
    end
    @testset "Inner transducer terminated during complete" begin
        @testset for xs in iterator_variants(1:5)
            @test xs |> xf |> Take(3) |> collect == [[1, 2], [3], [4, 5]]
        end
    end
end

@testset "Scan" begin
    @testset for xs in iterator_variants(1:10)
        xs isa Base.Generator && continue
        @test collect(Scan(+), xs) == cumsum(xs)
        @test collect(Scan((a, b) -> a + b), xs) == cumsum(xs)
    end

    xs0 = [0, -1, 3, -2, 1]
    @testset for xs in [xs0, collect(xs0)]
    # @testset for xs in iterator_variants([0, -1, 3, -2, 1])  # TODO: fix
        @test collect(Scan(max), xs) == [0, 0, 3, 3, 3]
        @test collect(Scan(min), xs) == [0, -1, -1, -2, -2]
    end
end

@testset "Scanx" begin
    @testset for xs in iterator_variants(1:10)
        xs isa Base.Generator && continue
        @test collect(Scanx(+), xs) == [0;cumsum(xs)]
        @test collect(Scanx(*), xs) == [1;cumprod(xs)]
        @test collect(Scanx((a, b) -> a + b), xs) == [0;cumsum(xs)]
    end

    xs0 = [0, -1, 3, -2, 1]
    @testset for xs in [xs0, collect(xs0)]
        # TODO: confirm that -Inf/Inf are the right monoidal identities here
        @test collect(Scanx(max), xs) == [-Inf,0, 0, 3, 3, 3]
        @test collect(Scanx(min), xs) == [Inf,0, -1, -1, -2, -2]
    end
end

@testset "ScanEmit" begin
    @testset for xs in iterator_variants(1:3)
        @test collect(ScanEmit(tuple, 0), xs) == 0:2
        @test collect(ScanEmit(tuple, 0, identity), xs) == 0:3
        @test collect(ScanEmit(tuple, nothing), xs) == [nothing; 1:2]
    end
    @testset "Inner transducer terminated during complete" begin
        @testset for xs in iterator_variants(1:3)
            @test xs |> ScanEmit(tuple, 0, identity) |> Take(4) |> collect == 0:3
        end
    end

    @testset "Do not call `complete` when reduced" begin
        xs = 1:8
        xf = opcompose(ScanEmit(CopyInit([]), identity) do u, x
            push!(u, x)
            if x % 3 == 0
                return u, []
            else
                return nothing, u
            end
        end, NotA(Nothing))

        @testset "foreach" begin
            called_with = []
            @test foreach(xf, xs) do chunk
                push!(called_with, copy(chunk))
                5 ∈ chunk && reduced(true)
            end == reduced(true)
            @test called_with == [1:3, 4:6]
        end

        @testset "foldl" begin
            called_with = []
            history = []
            @test foldl(xf, xs; init=false) do state, chunk
                push!(history, state)
                push!(called_with, copy(chunk))
                5 ∈ chunk && reduced(true)
            end == true
            @test called_with == [1:3, 4:6]
            @test history == [false, false]
        end
    end
end

@testset "ZipSource" begin
    @testset for xs in iterator_variants(1:5)
        @test collect(ZipSource(opcompose(Filter(isodd), Map(inc))), xs) ==
            collect(zip(1:2:5, 2:2:6))
    end

    @testset for xs in iterator_variants(1:6)
        @test xs |> Map(inc) |> ZipSource(Filter(isodd)) |> Map(first) |> collect == 3:2:7
    end

    @testset "Combination with stateful transducers" begin
        @testset for xs in iterator_variants(2:2:6)
            @test xs |> ZipSource(Map(identity)) |> Count() |> collect == 1:3
            @test xs |> ZipSource(Count()) |> Count() |> collect == 1:3
            @test xs |> Count() |> ZipSource(Map(x -> x + 10)) |> collect ==
                collect(zip(1:3, (1:3) .+ 10))
            @test (
                xs |>
                Enumerate() |>
                ZipSource(Map(x -> x[end] + 10)) |>
                Enumerate() |>
                collect
            ) == [
                (1, ((1, 2), 12))
                (2, ((2, 4), 14))
                (3, ((3, 6), 16))
            ]
        end
    end

    @testset "Union type input" begin
        @testset for xs in iterator_variants([10, missing, 20])
            @test collect(Zip(Count(), NotA(Missing)), xs) == [
                (1, 10)
                (3, 20)
            ]
        end
    end

    @testset "Nested ZipSource" begin
        @testset for xs in iterator_variants(1:5)
            @test collect(
                ZipSource(opcompose(Filter(isodd), Map(inc), ZipSource(Map(inc)))),
                xs,
            ) == [
                (1, (2, 3))
                (3, (4, 5))
                (5, (6, 7))
            ]
            @test collect(ZipSource(ZipSource(Map(inc))), xs) isa Vector
        end
    end
end

@testset "Replace" begin
    @testset for xs in iterator_variants(1:3)
        @test collect(Replace(Dict(1 => 10, 2 => 20)), xs) == [10, 20, 3]
    end
end

@testset "Take" begin
    @testset for xs in iterator_variants(1:5)
        @test xs |> Take(3) |> collect == 1:3
        @test xs |> Take(3) |> Take(4) |> collect == 1:3
        @test xs |> Take(4) |> Take(3) |> collect == 1:3
        @test xs |> Take(10) |> collect == 1:5
        @test xs |> Take(10) |> Take(20) |> collect == 1:5
        @test xs |> Take(20) |> Take(10) |> collect == 1:5
    end
    @testset "Combination with stateless transducers" begin
        @testset for xs in iterator_variants(1:5)
            @test xs |> Filter(isodd) |> Take(2) |> collect == 1:2:3
            @test xs |> Filter(iseven) |> Take(3) |> collect == 2:2:4  # no abort
            @test xs |> Take(3) |> Filter(isodd) |> collect == 1:2:3
            @test xs |> Take(3) |> Filter(iseven) |> collect == [2]
        end
    end
    @testset "Combination with stateful transducers" begin
        @testset for xs in Any[1:5, collect(1:5)]
            @test eltype(xs) <: Number
            @test xs |> Take(3) |> Scan(+) |> collect == cumsum(1:3)
            @test xs |> Scan(+) |> Take(3) |> collect == cumsum(1:3)
        end
        @testset for xs in iterator_variants(1:5)
            @test xs |> TakeLast(4) |> Take(2) |> collect == 2:3
            @test xs |> Take(4) |> TakeLast(2) |> collect == 3:4
        end
    end
end

@testset "TakeLast" begin
    @testset for xs in iterator_variants(1:5)
        @test xs |> TakeLast(2) |> collect == 4:5
        @test xs |> TakeLast(10) |> collect == 1:5
        @test xs |> TakeLast(3) |> TakeLast(4) |> collect == 3:5
        @test xs |> TakeLast(4) |> TakeLast(3) |> collect == 3:5
    end
    @testset "Combination with stateful transducers" begin
        @testset for xs in iterator_variants(1:5)
            @test xs |> Take(4) |> TakeLast(2) |> collect == 3:4
            @test xs |> TakeLast(4) |> Take(2) |> collect == 2:3
        end
    end
end

@testset "TakeWhile" begin
    @testset for xs in iterator_variants(1:5)
        @test xs |> TakeWhile(x -> x < 3) |> collect == 1:2
    end
    @testset "Combination with stateful transducers" begin
        @testset for xs in iterator_variants(1:5)
            @test xs |> TakeWhile(x -> x ≤ 4) |> TakeLast(2) |> collect == 3:4
            @test xs |> TakeLast(4) |> TakeWhile(x -> x ≤ 3) |> collect == 2:3
        end
    end
end

@testset "TakeNth" begin
    @testset for xs in iterator_variants(1:10)
        @test xs |> TakeNth(3) |> collect == 1:3:10
    end
    @testset "Combination with stateful transducers" begin
        @testset for xs in iterator_variants(1:10)
            @test xs |> Take(4) |> TakeNth(2) |> collect == (1:4)[1:2:end]
            @test xs |> TakeNth(2) |> Take(4) |> collect == (1:2:10)[1:4]
        end
    end
end

@testset "Drop" begin
    @testset for xs in iterator_variants(1:10)
        @test xs |> Drop(3) |> collect == 4:10
        @test xs |> Drop(3) |> Drop(3) |> collect == 7:10
    end
    @testset "Combination with stateful transducers" begin
        @testset for xs in iterator_variants(1:10)
            @test xs |> Drop(3) |> Take(3) |> collect == 4:6
            @test xs |> Take(5) |> Drop(3) |> collect == 4:5
        end
    end
end

@testset "DropLast" begin
    @testset for xs in iterator_variants(1:10)
        @test xs |> DropLast(3) |> collect == 1:7
        @test xs |> DropLast(3) |> DropLast(3) |> collect == 1:4
    end
    @testset "Combination with stateful transducers" begin
        @testset for xs in iterator_variants(1:10)
            @test xs |> Take(5) |> DropLast(3) |> collect == 1:2
            @test xs |> DropLast(3) |> Take(5) |> collect == 1:5
            @test xs |> DropLast(3) |> Take(10) |> collect == 1:7
        end
    end
end

@testset "DropWhile" begin
    @testset for xs in iterator_variants(1:10)
        @test xs |> DropWhile(x -> false) |> collect == 1:10
        @test xs |> DropWhile(x -> true) |> collect == []
        @test xs |> DropWhile(x -> x % 3 != 0) |> collect == 3:10
    end
    @testset "Combination with stateful transducers" begin
        @testset for xs in iterator_variants(1:10)
            @test xs |>
                  DropWhile(x -> x ÷ 3 == 0) |>
                  DropWhile(x -> x ÷ 3 == 1) |>
                  collect == 6:10
            @test xs |> DropWhile(x -> x ÷ 3 == 0) |> Take(3) |> collect == 3:5
            @test xs |> Take(3) |> DropWhile(x -> x ÷ 3 == 0) |> collect == [3]
            @test xs |> Take(2) |> DropWhile(x -> x ÷ 3 == 0) |> collect == []
        end
    end
end

@testset "FlagFirst" begin
    one2three = [
        (true, 1),
        (false, 2),
        (false, 3),
    ]
    @testset for xs in iterator_variants(1:3)
        @test xs |> FlagFirst() |> collect == one2three
    end
    @testset "Combination with stateful transducers" begin
        @testset for xs in iterator_variants(1:3)
            @test xs |> FlagFirst() |> FlagFirst() |> collect == [
                (true, (true, 1)),
                (false, (false, 2)),
                (false, (false, 3)),
            ]
            @test xs |> FlagFirst() |> Take(5) |> collect == one2three
            @test xs |> Take(5) |> FlagFirst() |> collect == one2three
            @test xs |> FlagFirst() |> Take(2) |> collect == one2three[1:2]
            @test xs |> Take(2) |> FlagFirst() |> collect == one2three[1:2]
        end
    end
end

# https://clojuredocs.org/clojure.core/keep
@testset "Keep" begin
    xf = @test_deprecated Keep() do x
        if x % 3 == 0
            x
        end
    end
    DEPWARN_ERROR || @testset for xs in iterator_variants(1:10)
        @test xs |> xf |> collect == 3:3:10
    end
end

@testset "KeepSomething" begin
    @test [Some(1), nothing, 2, Some(nothing)] |> KeepSomething() |> collect ==
          [1, 2, nothing]

    xf = KeepSomething() do x
        if x == 0
            :zero
        elseif x == 1
            Some(:one)
        end
    end

    @test collect(xf, 0:3) == [:zero, :one]
end

# https://clojuredocs.org/clojure.core/distinct
@testset "Unique" begin
    @testset for xs in iterator_variants([1, 1, 2, 1, 3, 2])
        @test xs |> Unique() |> collect == [1, 2, 3]
    end

    @testset for xs in iterator_variants([1, 1, 2, -1, 3, 2])
        @test xs |> Unique(x -> x^2) |> collect == [1, 2, 3]
    end
end

# https://clojuredocs.org/clojure.core/interpose
@testset "Interpose" begin
    @testset for xs in iterator_variants(1:3)
        @test xs |> Interpose(0) |> collect == [1, 0, 2, 0, 3]
    end
end

# https://clojuredocs.org/clojure.core/dedupe
@testset "Dedupe" begin
    @testset for xs in iterator_variants([1, 1, 1, 2, 3, 3, 1])
        @test xs |> Dedupe() |> collect == [1:3; 1]
    end
end

# https://clojuredocs.org/clojure.core/partition-all
@testset "Partition" begin
    @testset for xs in iterator_variants(1:10)
        @test xs |> Partition(3) |> Map(copy) |> collect == [[1:3;], [4:6;], [7:9;]]
    end
    @testset for xs in iterator_variants(1:6)
        @test xs |> Partition(3, 1) |> Map(copy) |> collect ==
            [[1:3;], [2:4;], [3:5;], [4:6;]]
    end
    @testset "w=3, 7:$i" for i in 7:9
        @testset for xs in iterator_variants(1:i)
            @test xs |> Partition(3, flush=true) |> Map(copy) |> collect ==
                [[1:3;], [4:6;], [7:i;]]
        end
    end
    @testset "w=4, 5:$i" for i in 5:8
        @testset for xs in iterator_variants(1:i)
            @test xs |> Partition(4, flush=true) |> Map(copy) |> collect ==
                [[1:4;], [5:i;]]
        end
    end
    @testset "w=4, 9:$i" for i in 9:12
        @testset for xs in iterator_variants(1:i)
            @test xs |> Partition(4, flush=true) |> Map(copy) |> collect ==
                [[1:4;], [5:8;], [9:i;]]
        end
    end
    @testset "Combination with stateful transducers" begin
        @testset for xs in iterator_variants(1:6)
            @test xs |> Take(5) |> Partition(3, 1) |> Map(copy) |> collect ==
                [[1:3;], [2:4;], [3:5;]]
            @test xs |> Partition(3, 1) |> Map(copy) |> Take(3) |> collect ==
                [[1:3;], [2:4;], [3:5;]]
        end
    end
    @testset "Inner transducer terminated during complete" begin
        @testset "w=4, 9:$i" for i in 9:12
            @testset for xs in iterator_variants(1:i)
                @test xs |> Partition(4, flush = true) |> Map(copy) |> Take(3) |> collect ==
                      [[1:4;], [5:8;], [9:i;]]
            end
        end
    end
end

@testset "Iterated" begin
    @testset for xs in iterator_variants(1:3)
        @test xs |> Iterated(x -> x + 1, 1) |> collect == 1:3
    end
    @testset "Combination with stateful transducers" begin
        @testset for xs in iterator_variants(1:3)
            @test xs |> Iterated(identity, nothing) |> Iterated(x -> x + 1, 1) |> collect ==
                  1:3
            @test xs |> Iterated(x -> x + 1, 1) |> Drop(1) |> collect == 2:3
            @test xs |> Drop(1) |> Iterated(x -> x + 1, 1) |> collect == 1:2
        end
    end
end

@testset "Count" begin
    @testset for xs in iterator_variants(1:3)
        @test xs |> Count() |> collect == 1:3
        @test xs |> Count(1im) |> collect == [0+1im, 1+1im, 2+1im]
        @test xs |> Count(Second(1)) |> collect == Second.(xs)
        @test_broken xs |> Count(Second(1), Millisecond(1)) |> collect == [
            Second(1) + Millisecond(0),
            Second(1) + Millisecond(1),
            Second(1) + Millisecond(2),
        ]
    end
    @testset "Combination with stateful transducers" begin
        @testset for xs in iterator_variants(1:3)
            @test xs |> Count() |> Drop(1) |> collect == 2:3
            @test xs |> Drop(1) |> Count() |> collect == 1:2
        end
    end
end

@testset "(Get|Set)Index" begin
    for T in [GetIndex, SetIndex]
        @test T([0 + 1im]) == T([0 + 1im])
        @test T{true}([0]) != T([0])
        @test T([0]) != T{true}([0])
        @test T([0]) != T([0im])  # should it?
    end
end

@testset "Inject" begin
    @testset for xs in iterator_variants(1:3)
        @test xs |> Inject(xs) |> collect == collect(zip(xs, xs))
        @test xs |> Inject(5:6) |> collect == collect(zip(1:2, 5:6))
        @test 5:6 |> Inject(xs) |> collect == collect(zip(5:6, 1:2))
    end
    @testset "Combination with stateful transducers" begin
        @testset for xs in iterator_variants(1:3)
            @test xs |> Inject(xs) |> Take(2) |> collect == collect(zip(1:2, 1:2))
            @test xs |> Take(2) |> Inject(xs) |> collect == collect(zip(1:2, 1:2))
            @test xs |> Inject(xs) |> TakeLast(2) |> collect == collect(zip(2:3, 2:3))
            @test xs |> TakeLast(2) |> Inject(xs) |> collect == collect(zip(2:3, 1:2))
            @test xs |> Inject(1:1) |> TakeLast(2) |> collect == collect(zip(1:1, 1:1))
            @test xs |> TakeLast(2) |> Inject(1:1) |> collect == collect(zip(2:2, 1:1))
        end
    end
end

@testset "Enumerate" begin
    @testset for xs in iterator_variants(2:2:6)
        @test xs |> Enumerate() |> collect == [(1,2), (2,4), (3,6)]
        @test xs |> Enumerate(10, 20) |> collect == [(10, 2), (30, 4), (50, 6)]
    end
    @testset "Combination with stateful transducers" begin
        @testset for xs in iterator_variants(2:2:6)
            @test xs |> Enumerate() |> Take(2) |> collect == [(1, 2), (2, 4)]
            @test xs |> Drop(2) |> Enumerate() |> collect == [(1, 6)]
        end
    end
end

@testset "OfType" begin
    @testset "_next_oftype for non-tuple" begin
        @testset for xs in iterator_variants([0.0, 1.0, missing])
            @test xs |> OfType(Float64) |> collect == [0.0, 1.0]
            @test xs |> OfType(Number) |> collect == [0.0, 1.0]
        end
        @testset for xs in iterator_variants(["a", :b])
            @test xs |> OfType(Symbol) |> collect == [:b]
            @test xs |> OfType(AbstractString) |> collect == ["a"]
        end
    end
    @testset "_next_oftype for 2-tuple" begin
        @testset for xs in iterator_variants([(0.0, missing),
                                              (missing, 0.0),
                                              (1.0, 2.0)])
            @test xs |> OfType(Tuple{Float64,Float64}) |> collect == [(1.0, 2.0)]
            @test xs |> OfType(Tuple{Number,Float64}) |> collect == [(1.0, 2.0)]
            @test xs |> OfType(Tuple{Float64,Number}) |> collect == [(1.0, 2.0)]
            @test xs |> OfType(Tuple{Vararg{Float64}}) |> collect == [(1.0, 2.0)]
            @test xs |> OfType(Tuple{Vararg{Number}}) |> collect == [(1.0, 2.0)]
        end
    end
    @testset "_next_oftype_t" begin
        @testset "$(typeof(xs))" for xs in iterator_variants([
                [Tuple(i == j ? missing : i for i in 1:10) for j in 1:10]
                Tuple(1:10)
                ])
            desired = [Tuple(1:10)]
            @test xs |> OfType(Tuple{Vararg{Int}}) |> collect == desired
            @test xs |> OfType(Tuple{Vararg{Number}}) |> collect == desired
        end
    end
end

@testset "ReduceIf" begin
    @test foldl(right, ReduceIf(x -> x == 3), 1:10) === 3
    @test foldl(right, ReduceIf(x -> x == 3), 1:2) === 2
end

@testset "AbortIf" begin
    @test foldl(right, AbortIf(x -> x == 3), 1:10) === 2
    @test foldl(right, AbortIf(x -> x == 3), 1:1) === 1
end

@testset "Invalid arguments" begin
    @test_throws ArgumentError Take(-1)
    @test_throws ArgumentError TakeLast(-1)
    @test_throws ArgumentError TakeNth(0)
    @test_throws ArgumentError Drop(-1)
    @test_throws ArgumentError DropLast(-1)
    @test_throws ArgumentError Partition(0)
    @test_throws ArgumentError Partition(0, 1)
    @test_throws ArgumentError Partition(1, 0)
    @test_throws ArgumentError TCat(0)
end

end  # module