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

using MXNet, Plots, StatsBase, LaTeXStrings

#pyplot() # workaround for Julia 0.6 Plots.jl bug

export saveplot
function saveplot(env::DLEnv, name; formats=["png", "pdf", "svg"])
  for format in formats
    savefig(joinpath(env.dir, "plots", name*"."*format))
  end
  return plot!()
end

export log_histogram2d
function log_histogram2d(x, y, x_bins, y_bins; colorbar=:none)
  mse_hist_2d = fit(Histogram{Float64}, (x, y), (x_bins, y_bins), closed=:left)
  broadcast!(x -> x <= 0 ? NaN : log10(x), mse_hist_2d.weights, mse_hist_2d.weights)
  plot(mse_hist_2d, colorbar=colorbar, framestyle=:box, color=:viridis)
end

export paperplot
function paperplot(ratio=1.2, columns=1)
    gr()
    f = 10
    width = 86 * columns * 4.2524
    height = width / ratio
    plot(size=(width, height), dpi=66.666, xtickfontsize=f, ytickfontsize=f, legendfontsize=f, xguidefontsize=f, yguidefontsize=f, titlefontsize=f, framestyle=:box)
end

export plot_waveforms
"""
Plots a number of waveforms from the given EventLibrary in a single diagram. The figure is saved in the given environment using the name of the library.
"""
function plot_waveforms(env::DLEnv, events::EventLibrary; count=4, bin_width_ns=10, cut=nothing, diagram_font=font(16), evt_indices=nothing)
  filepath = joinpath(env.dir, "plots", "waveforms-$(events[:name]).png")
  info(env, 2, "Plotting waveforms to $filepath")

  if evt_indices == nothing
    count = min(count, eventcount(events))
    data = waveforms(events)[:,1:count]
  else
    data = waveforms(events)[:, evt_indices]
  end

  if cut != nothing
    data = data[cut,:]
  end

  x_axis = linspace(0, (size(data,1)-1)*bin_width_ns, size(data,1))
  plot(size=(600, 400), legend=:none)
  for i in 1:size(data, 2)
    plot!(x_axis, data[:,i], linewidth=3)
  end
  xaxis!("Time (ns)", diagram_font)
  yaxis!("Current pulse", diagram_font)
  savefig(filepath)

  return events
end

"""
For each EventLibrary, a number of waveforms are plotted in one diagram. All figures are saved in the given environment using the names of the libraries.
"""
function plot_waveforms(env::DLEnv, data::DLData; count=4, bin_width_ns=10, cut=nothing)
  for lib in data
    plot_waveforms(env, lib; count=count, bin_width_ns=bin_width_ns, cut=cut)
  end
  return data
end

export plot_waveform_thumbnail
function plot_waveform_thumbnail(env::DLEnv, waveform::Vector{Float32}; bin_width_ns=10, tcut=50:200, titlestring=nothing, diagram_size=(140, 100), diagram_font=font(12), filename="waveform-thumbnail.png")
  waveform = waveform[tcut]
  tlength = tcut[end] - tcut[1]
  time = linspace(0, (tlength-1)*0.01, length(waveform))

  fig = plot(time, waveform, size=diagram_size, line=(:blue), legend=:none, xformatter = x -> "$(Int(x)) us", grid=false)
  yticks!(Float64[])
  xticks!([0,1])
  xaxis!(diagram_font)
    valuerange = maximum(waveform) - minimum(waveform)
  yaxis!((minimum(waveform)-valuerange*0.1, maximum(waveform)+valuerange*0.1))
  if titlestring != nothing
    title!(titlestring, titlefont=diagram_font)
  end
  savefig(joinpath(env.dir, "plots", filename))
  fig
end

export waveform_thumbnail
function waveform_thumbnail(events::EventCollection, waveform, color; tcut=50:200, diagram_font=font(10), diagram_size_mm=(40,40), energy=nothing)
    waveform = waveform[tcut]
    tlength = tcut[end] - tcut[1]
    time = linspace(0, (tlength-1)*sampling_period(events)*1e6, length(waveform))

    if energy != nothing
        waveform *= energy / maximum(waveform)
        diagram_size_mm = (diagram_size_mm[1] * 1.25, diagram_size_mm[2])
    end

    fig = plot(time, waveform, size=diagram_size_mm .* 4.2524, dpi=66.666, line=color, legend=:none, xformatter = x -> "$((x==0)?0:x) \\mus", grid=false, framestyle=:box)

    xticks!(Float64[])
    xaxis!(diagram_font)

    valuerange = maximum(waveform) - minimum(waveform)

    yaxis!((-valuerange*0.1, maximum(waveform)+valuerange*0.1), diagram_font)
    if energy == nothing
        yticks!(Float64[])
    else
        yaxis!("Energy (keV)")
    end
    fig
end

# NETWORK LEARNING

function plot_learning_curves(n::NetworkInfo, filename; from_zero::Bool=false)
  train = n.training_curve
  xval = n.xval_curve
  plot(size=(800,500), legendfont=font(20))
  plot!(train, label="Training set", line=2)
  plot!(xval, label="Validation set", line=2)
  xaxis!("Training time (epochs)", font(20))
  yaxis!("Mean squared error", font(20))
  if from_zero
    yaxis!((0, max(maximum(train), maximum(xval))))
  end
  savefig("$filename.pdf")
  savefig("$filename.png")
end

export plot_dnn
function plot_dnn(env::DLEnv, n::NetworkInfo)
  dir = joinpath(env.dir, "plots", n.name)
  isdir(dir) || mkdir(dir)
  info(env, 2, "Generating plots in $dir...")
  plot_learning_curves(n, joinpath(dir, "learning_curves"))
end

# AUTOENCODER

export plot_autoencoder
function plot_autoencoder(env::DLEnv, n::NetworkInfo)
  dir = joinpath(env.dir, "plots", n.name)
  isdir(dir) || mkdir(dir)
  info(env, 2, "Generating plots in $dir...")

  model = n.model

  plot_learning_curves(n, joinpath(dir, "learning_curves"); from_zero=true)
  visualize_1D_convolution(model, :conv_1_weight, joinpath(dir,"filters1.png"))
  visualize_2D_convolution(model, :conv_2_weight, joinpath(dir,"filter2"))
  info(env, 2, "Saved network plots to $dir")
  return model, dir
end

function plot_autoencoder(env::DLEnv, n::NetworkInfo, data::Dict{Symbol,EventLibrary}; count=20, transform=identity)
  for (key,value) in data
    plot_autoencoder(env, n, value; count=count, transform=transform)
  end
end

function plot_autoencoder(env::DLEnv, n::NetworkInfo, data::EventCollection; count=20, transform=identity)
  batch_size = n["batch_size"]

  model, dir = plot_autoencoder(env, n)

  for lib in data

      if eventcount(lib) > count
        lib = lib[1:count]
      end

      provider = padded_array_provider(:data, waveforms(lib), batch_size)
      plot_reconstructions([n.model], ["Reconstruction"], waveforms(lib), provider, dir, file_prefix=lib[:name]*"-", sample_count=eventcount(lib), transform=transform, bin_width_ns=sampling_period(lib), legendpos=(lib[:waveform_type] == "current" ? :topright : :topleft))
    end

  # loss = model.arch
  # open(joinpath(dir,"graphviz.dot"), "w") do file
  #   info(file, mx.to_graphviz(loss))
  # end

  info(env, 2, "Saved reconstruction plots to $dir")

  return model, dir
end

function plot_reconstructions(models, names, waveforms, provider, plot_dir;
    file_prefix="", sample_count=20, transform=identity,
    bin_width_ns=10, y_label="Current pulse", diagram_font=font(14), legendpos=:topright
    )
  # plot_waveform_comparisons(encode_decode(...), ...)
  pred_evals = [transform(mx.predict(model, provider)) for model in models]
  truth = transform(waveforms)

  sample_count = min(sample_count, size(truth, 2))
  x_axis = linspace(0, (size(truth,1)-1)*bin_width_ns*1e9, size(truth,1))

  reconst_animation = @animate for i in 1:sample_count
    plot(x_axis, truth[:,i], label="Data", linewidth=2, legendfont=diagram_font, legend=legendpos)
    for modelI in 1:length(models)
      plot!(x_axis, pred_evals[modelI][:,i], linewidth=3, label=names[modelI])
    end
   xaxis!("Time (ns)", diagram_font)
   yaxis!(y_label, diagram_font)
   savefig("$plot_dir/reconst-$(file_prefix)$i.png")
  end
  gif(reconst_animation, "$plot_dir/reconstructions.gif", fps=1)
end

export plot_waveform_comparisons
function plot_waveform_comparisons(env::DLEnv, libs::EventLibrary...; count=20, cut=nothing, diagram_font=font(14), transform=identity, title::String="", y_label="Current pulse")
  diag_count = min(count, minimum(length.(libs)))
  if diag_count != count
    info(env, 2, "Diagram count reduced to $diag_count, requested: $count")
  end
  sample_count = minimum(samples.(libs))
  x_axis = linspace(0, (sample_count-1)*sampling_period(libs[1])*1e9, sample_count)

  # Create directory
  dir = joinpath(env, "plots", title)
  isdir(dir) || mkdir(dir)
  info(env, 2, "Saving comparison plots to $dir ($(length(libs)) libraries)")

  for i in 1:diag_count
    plot(legendfont=diagram_font, legend=:topright)
    for lib in libs
      plot!(x_axis, waveforms(lib)[:,i], linewidth=2, label=lib[:name])
    end
   xaxis!("Time (ns)", diagram_font)
   yaxis!(y_label, diagram_font)

   savefig(joinpath(dir, "wfcmp-$(join(name.(libs),"-"))$i.png"))
  end
  return libs
end

export visualize_1D_convolution
function visualize_1D_convolution(model, name, filename)
  pyplot_config()

  colors = [:blue :red :green :purple]
  for param in model.arg_params
    if(param[1] == name)
      filters = copy(param[2])[1,:,1,:] # One column per filter
      filters = vcat(filters, transpose(filters[end,:]))
      plot(filters, line=(2, :steppost), label="", color=colors)
      title!("Learned convolutional filters in first layer")
      xaxis!("Sample")
      yaxis!("Conv. filter value")
      savefig(filename)
    end
  end
end

function visualize_2D_convolution(model, name, filename_base)
  for param in model.arg_params
    if(param[1] == name)
      for i in 1:size(param[2], 4)
        filters = copy(param[2])[1,:,:,i]
        heatmap(filters)
        title!("Learned convolutional filters in second layer")
        xaxis!("Delta time")
        yaxis!("Feature map")
        savefig("$filename_base-$i.png")
      end
    end
  end
end

export plot_reconstruction_accuracy
function plot_reconstruction_accuracy(env, net::NetworkInfo, pulses::EventCollection, noise::EventCollection; waveform_scale::Real=256, filename=nothing, include_examples=false, reconst=nothing, E_axis=linspace(500, 2700, 250), err_axis=nothing)
  if reconst == nothing
    info(env, 3, "Computing pulse reconstructions...")
    reconst = scale_waveforms(encode_decode(scale_waveforms(pulses, waveform_scale), net), 1.0/waveform_scale)
  end

    info(env, 3, "Calculating standard deviations...")
    std_reconst = zeros(Float32, eventcount(noise))
    std_noise = zeros(Float32, eventcount(noise))

    wf_pulses = waveforms(pulses)
    wf_reconst = waveforms(reconst)
    wf_noise = waveforms(noise)

    for i in 1:eventcount(pulses)
        std_reconst[i] = std(wf_reconst[:,i]-wf_pulses[:,i])
        std_noise[i] = std(wf_noise[:,i])
    end

    avg_std = std(std_reconst-std_noise)
    info(env, 3, "Std of additional reconstruction deviation: $avg_std")

    # Plot added reconstruction error
    if err_axis == nothing
      err_axis = linspace(-3avg_std+0.0005, 0.0015+2avg_std, 200)
    end
    info(env, 3, "Fitting 2D histogram with $(length(std_reconst)) entries...")
    mse_hist_2d = fit(Histogram{Float64}, (pulses[:E], std_reconst-std_noise), (E_axis, err_axis), closed=:left)
    broadcast!(x -> x <= 0 ? NaN : log10(x), mse_hist_2d.weights, mse_hist_2d.weights)
    info(env, 3, "Plotting 2D histogram...")
    fig_err = plot(mse_hist_2d, colorbar=:none, framestyle=:box)
    xaxis!("Energy (keV)", xticks=[500, 1000, 1500, 2000, 2500])
    yaxis!(L"$\sigma_{reconst} - \sigma_{noise}$")

    info(env, 2, "Average added error in terms of charge: $(median(std_reconst-std_noise))")
    #println("Average added error in terms of A/E: ", median((std_reconst-std_noise) ./ AoE))
    #fig_err_hist = stephist((std_reconst-std_noise), bins=y_axis)
    mse_hist_1d = fit(Histogram{Float64}, std_reconst-std_noise, err_axis, closed=:left)
    he_indices = find(E->E>1000, pulses[:E])
    info(env, 3, "Fitting 1D histograms...")
    mse_hist_1d_he = fit(Histogram{Float64}, std_reconst[he_indices]-std_noise[he_indices], err_axis, closed=:left)
    fig_err_hist = plot(mse_hist_1d.weights/1000, err_axis, line=(:black, :steppre), label="All events", framestyle=:box)
    plot!(mse_hist_1d_he.weights/1000, err_axis, line=(:lightblue, :steppre), label="> 1 MeV")
    xaxis!("Events (1000)")
    yaxis!(yticks=nothing)

    info(env, 2, "Plotting results to $(joinpath(env, "plots", net.name, "Added noise.pdf"))")
    plot(fig_err, fig_err_hist, size=(800, 400), layout=@layout([a{0.75w} b]))
    if filename == nothing
        savefig(joinpath(env, "plots", net.name, "Added noise.pdf"))
        savefig(joinpath(env, "plots", net.name, "Added noise.png"))
    end
    if filename != nothing && !include_examples
        savefig(filename)
    end

    if include_examples
      # Plot examples of SSE and MSE
      SSE = SSE_at(pulses, 2038) # 1592.5
      MSE = MSE_at(pulses, 2038) # 1620.7
      info(env, 3, "Drawing SSE at $(pulses[:E][SSE]) and MSE at $(pulses[:E][MSE])")

      fig_SSE = reconstruction_plot(pulses, reconst, SSE, "Single-site event", true; zoom=50:206)
      fig_MSE = reconstruction_plot(pulses, reconst, MSE, "Multi-site event", false; zoom=50:206)

      all_layout = @layout [grid(1,2,widths=[0.50, 0.5]){0.45h}
                          a{0.7w} b]
      plot(fig_SSE, fig_MSE, fig_err, fig_err_hist, size=(15 * 39.37, 15 * 39.37), layout=all_layout)

      if filename == nothing
          savefig(joinpath(env, "plots", net.name, "reconstructions_$(pulses[:detector_name]).pdf"))
          savefig(joinpath(env, "plots", net.name, "reconstructions_$(pulses[:detector_name]).png"))
      end

      if filename != nothing
          savefig(filename)
      end
    end

    return pulses, noise, reconst
end

export reconstruction_accuracy_percent_fig
function reconstruction_accuracy_percent_fig(pulses::EventCollection, reconst::EventCollection, noise::EventCollection; E_axis=500:10:2700, err_axis=0.75:0.01:1.8)
    std_reconst, std_noise = calculate_deviation(pulses, reconst, noise)

    # Plot added reconstruction error
    if err_axis == nothing
      err_axis = linspace(0.75, 2, 100)
    end
    mse_hist_2d = fit(Histogram{Float64}, (pulses[:E], std_reconst./std_noise), (E_axis, err_axis), closed=:left)
    broadcast!(x -> x <= 0 ? NaN : log10(x), mse_hist_2d.weights, mse_hist_2d.weights)
    fig_err = plot(mse_hist_2d, colorbar=:none, framestyle=:box, color=:viridis)
    xaxis!("Energy (keV)", xticks=[500, 1000, 1500, 2000, 2500])
    yaxis!(L"$\sigma_{reconst} / \sigma_{noise}$")
    plot!(x->1, line=(:white, 0.5, :dash), label="")

    mse_hist_1d = fit(Histogram{Float64}, std_reconst./std_noise, err_axis, closed=:left)
    he_indices = find(E->E>1000, pulses[:E])
    mse_hist_1d_he = fit(Histogram{Float64}, std_reconst[he_indices]-std_noise[he_indices], err_axis, closed=:left)
    fig_err_hist = plot(mse_hist_1d.weights/1000, err_axis[2:end], line=(:black, :steppre), label="All events", framestyle=:box)
    plot!(mse_hist_1d_he.weights/1000, err_axis[2:end], line=(:lightblue, :steppre), label="> 1 MeV")
    xaxis!("Events (1000)")
    yaxis!(yticks=nothing)
    plot!(x->1, line=(:black, 0.5, :dash), label="")

    plot(fig_err, fig_err_hist, size=(800, 400), layout=@layout([a{0.75w} b]))
end

export reconstruction_plot
function reconstruction_plot(pulses::EventCollection, reconst::EventCollection, index::Integer, title::AbstractString, plot_y_axis; plot_x_axis=true, zoom=:)
    t_axis = collect(sample_times(pulses) * 1e6)[zoom]
    t_axis -= t_axis[1]
    fig = plot(t_axis, waveforms(pulses)[zoom,index], line=(1,:black), label="Data", framestyle=:box)
    plot!(t_axis, waveforms(reconst)[zoom,index], line=(2,:blue), label="Reconstruction")
    if plot_x_axis
        xaxis!(L"Time ($\mu$s)")
    end
    if plot_y_axis
        yaxis!("Current pulse")
    end
    yaxis!((min(minimum(waveforms(pulses)[zoom,index]), minimum(waveforms(reconst)[zoom,index])), 1.25 * max(maximum(waveforms(pulses)[zoom,index]), maximum(waveforms(reconst)[zoom,index]))))
    if length(title) > 0
      title!("$title ($(Int(round(pulses[:E][index]))) keV)", titlefont=font(11))
    end

    return fig
end

export energy_dependence_figure
function energy_dependence_figure(env::DLEnv, pulses::EventCollection, detector::AbstractString; cut_value="90% DEP")
    effs = get_peak_efficiencies(env, flatten(pulses), detector);
    cut_value = find_cut_value(effs, cut_value)
    info(env, 3, "Cut value: $cut_value")
    pass = filter(pulses, :psd, psd->psd>=cut_value)
    reject = filter(pulses, :psd, psd->psd<cut_value)
    info(env, 3, "Overall $(eventcount(pass)) passing, $(eventcount(reject)) rejected")

    equal_count_energy_edges = vcat(collect(equal_event_count_edges(pulses, :E)), 3000)

    energy_axis = equal_count_energy_edges# linspace(500, 3000, 200)
    hist_pass = fit(Histogram{Int64}, pass[:E], energy_axis, closed=:left)
    hist_reject = fit(Histogram{Int64}, reject[:E], energy_axis, closed=:left)
    hist_all = fit(Histogram{Float64}, pulses[:E], linspace(500, 3000, 250), closed=:left)
    fractions = hist_pass.weights ./ (hist_pass.weights.+hist_reject.weights)

    energy_lim = (500, 2650)
    fig_fraction = plot(framestyle=:box)
    plot!(hist_all.edges, hist_all.weights / maximum(hist_all.weights), line=(:steppost, :grey), fill=(0,:grey), label="Events (normalized)")
    plot!(hist_pass.edges, fractions, line=(:steppost, :blue), label="Passing fraction")
    xaxis!("Energy (keV)", energy_lim, xticks=[500, 1000, 1500, 2000, 2500])
    yaxis!("Fraction of events")

    e_dep_hist = fit(Histogram{Float64}, (convert(Array{Float64},pulses[:E]),
    convert(Array{Float64},pulses[:psd])), (linspace(500, 3000, 101), linspace(0, 1, 81)), closed=:left)
  broadcast!(x -> x <= 0 ? NaN : log10(x), e_dep_hist.weights, e_dep_hist.weights)
    fig_hist_2d = plot(e_dep_hist, framestyle=:box)
    xaxis!(energy_lim, xticks=nothing)
    yaxis!("Classifier output", colorbar=nothing)

    hist_class = fit(Histogram{Int64}, pulses[:psd], linspace(0, 1, 50), closed=:left)
    hist_class_ROI = fit(Histogram{Int64}, filter(pulses, :E, E->(E>1700)&&(E<2090))[:psd], linspace(0, 1, 100), closed=:left)
    fig_hist_1d = plot(yticks=nothing, xticks=nothing, framestyle=:box)
    plot!(hist_class.weights/maximum(hist_class.weights), hist_class.edges, line=(:steppre, :black), label="All events")
    plot!(hist_class_ROI.weights/maximum(hist_class_ROI.weights), hist_class_ROI.edges, line=(:steppre, :lightblue), label="ROI")

    empty_plot = plot(framestyle=:none)
    l = @layout([a{0.7w} b; c{0.4h} d])
    return plot(fig_hist_2d, fig_hist_1d, fig_fraction, empty_plot, layout=l, size=(15 * 39.37, 12 * 39.37))
end

# CLASSIFIER

export plot_classifier_histogram
function plot_classifier_histogram(dir, events::EventLibrary, label_key, psd_key;
    diagram_font=font(14), nbins = 30)

  if !haskey(events, psd_key)
    throw(ArgumentError("events must have label 'psd'."))
  end

  events = filter(events, :E, e -> (e < 3800 && e > 600))

  bins = linspace(minimum(events.labels[psd_key]), maximum(events.labels[psd_key]), nbins)

    # Labeled distribution
  if haskey(events, label_key)
    labels = events[label_key]
    psd = events[psd_key]

    psd_SSE = psd[find(l -> l==1, labels)]
    psd_MSE = psd[find(l -> l==0, labels)]

    # Overall PSD distribution
    histogram(psd_SSE, bins=bins, label="Tl DEP ($(length(psd_SSE)) events)", legendfont=diagram_font, linewidth=0)
    histogram!(psd_MSE, bins=bins, label="Bi FEP ($(length(psd_MSE)) events)", fillalpha=0.7, legendfont=diagram_font, linewidth=0)
    xaxis!("Classifier response", diagram_font)
    yaxis!("Event count", diagram_font)
    savefig(joinpath(dir,"Class_distribution_$(events[:name]).png"))
  end

    # Total distribution
  histogram(events[psd_key], bins=bins, label="All ($(eventcount(events)) events)", legendfont=diagram_font, linewidth=0)
  xaxis!("Classifier response", diagram_font)
  yaxis!("Event count", diagram_font)
  savefig(joinpath(dir,"Distribution $(events[:name]).png"))

    # Energy dependence
  e_dep_hist = fit(Histogram{Float64}, (convert(Array{Float64},events[:E]),
    convert(Array{Float64},events[psd_key])), (linspace(1000, 3000, 101), linspace(0, 1, 51)))
  broadcast!(x -> x <= 0 ? NaN : log10(x), e_dep_hist.weights, e_dep_hist.weights)
  plot(e_dep_hist, color=:viridis)
  xaxis!("Energy (keV)")
  yaxis!("Classifier output")
  title!("Classification vs energy (log10)")
  savefig(joinpath(dir, "Energy distributions $(events[:name]).png"))

    # Time dependence
    if haskey(events, :timestamp)
        events_sorted = sort(events, :timestamp)
      histogram2d(convert(Array{Float64}, 1:length(events_sorted[psd_key]))/1000, convert(Array{Float64},events_sorted[psd_key]))
        xaxis!("Entries (1000)", font(16))
        yaxis!("Classification", font(16))
        title!("Time stability of $(events[:detector_name])", titlefont=font(16))
      savefig(joinpath(dir, "Distribution over time $(events[:name]) entries.png"))
      histogram2d(convert(Array{Float64}, (events_sorted[:timestamp]-events_sorted[:timestamp][1])/60/60/24), convert(Array{Float64},events_sorted[psd_key]))
        xaxis!("Time (days)", font(16))
        yaxis!("Classification", font(16))
        title!("Time stability of $(events[:detector_name])", titlefont=font(16))
      savefig(joinpath(dir, "Distribution over time $(events[:name]) millis.png"))
    end

#  if haskey(events, :multiplicity)
#    histogram2d(convert(Array{Float64},events.labels[:multiplicity]), convert(Array{Float64},events.labels[psd_key]))
#    savefig(joinpath(dir, "Multiplicity correlation $(events[:name]).png"))
#  end

  return bins
end

export plot_energy_histogram
function plot_energy_histogram(env::DLEnv, events::EventLibrary, psd_cut, psd_key)
  all = copy(events)
  all.prop[:diagram_label] = "All events"
  cut = filter(all, psd_key, psd -> psd > psd_cut)
  cut.prop[:diagram_label] = "After PSD, cut value $(@sprintf("%.2f",psd_cut))"
  # range = linspace(1500, 2800, 150)
  plot_energy_histogram([all, cut], 150,
      joinpath(env.dir, "plots", "PSD-energies-cut$(@sprintf("%.3f",psd_cut)).png"))
end

function plot_energy_histogram(eventlibs::Array{EventLibrary}, nbins, filename;
  diagram_font=font(14), label=:diagram_label)

  histogram(xticks=[2000, 2500], legend=:bottomleft, legendfont=diagram_font)

  # Peaks before and after
  for events in eventlibs
    histogram!(events[:E], nbins=nbins, label=events[label], linewidth=0)
  end
  xaxis!("Energy (keV)", diagram_font)
  yaxis!("Events", :log10, (1,Inf), diagram_font)
  savefig(filename)
end

export efficiency_curve_figure
function efficiency_curve_figure(effs::EfficiencyCollection; invert_indices::Vector{Int}=Int[])
  cut_values = effs.curves[1].cut_values

  plot(size=(600,400), legendfont=font(16))
  for i in 1:length(effs.curves)
    curve = effs.curves[i]
    if i in invert_indices
      values = 1 .- curve.efficiencies
    else
      values = curve.efficiencies
    end
    plot!(curve.cut_values, values, line=(1), marker=(:circle), markerstrokewidth=0, label=replace(string(curve.name), "_", " "))
  end

  xaxis!("Cut value", font(20))
  yaxis!("Efficiency", (0,1), font(20))
end

export plot_efficiency_curves
function plot_efficiency_curves(file::AbstractString, effs::EfficiencyCollection)
  efficiency_curve_figure(effs)
  savefig("$file.png")
  savefig("$file.pdf")
  return effs
end

export plot_efficiencies
function plot_efficiencies(file::AbstractString, effs::EfficiencyCollection...)
  fig = plot(size=(700, 600), legendfont=font(16), legend=:bottomleft)

  for eff in effs
    if length(eff.curves[1].cut_values) <= 2
      plot!(roc_curve(eff), label=eff.title, line=(:dash, :grey, 2))
    else
      plot!(roc_curve(eff), label=eff.title, linewidth=4)
    end
  end
  xaxis!("Signal efficiency", (0.6,1), font(20))
  yaxis!("Background rejection", (0,1), font(20))
  savefig("$file.png")
  savefig("$file.pdf")
  return fig
end

export roc_plot
function roc_plot(effs::EfficiencyCollection...)
  fig = plot(size=(700, 600), legendfont=font(16), legend=:bottomleft)

  for eff in effs
    if length(eff.curves[1].cut_values) <= 2
      plot!(roc_curve(eff), label=eff.title, line=(:dash, :grey, 2))
    else
      plot!(roc_curve(eff), label=eff.title, linewidth=1, marker=:circle, markerstrokewidth=0)
    end
  end
  xaxis!("Signal efficiency", (0.6,1), font(20))
  yaxis!("Background rejection", (0,1), font(20))
  return fig
end

export plot_classifier
function plot_classifier(env::DLEnv, dirname::AbstractString, libs::EventLibrary...;
    classifier_key=:psd, label_key=:SSE, plot_AoE=false)
  dir = joinpath(env.dir, "plots", dirname)
  isdir(dir) || mkdir(dir)

  for lib in libs
    plot_classifier_histogram(dir, lib, label_key, classifier_key)
  end

  peak_effs = [get_peak_efficiencies(env, lib, lib[:detector_name], join(get_classifiers(lib), " + ")) for lib in libs]
  current_curve = load_current_effs(libs[1][:detector_name])

  for (i,lib) in enumerate(libs)
    plot_efficiency_curves(joinpath(dir, "PSD Efficiencies vs Cut $(lib[:name])"), peak_effs[i])
  end
  if current_curve != nothing
    plot_efficiency_curves(joinpath(dir, "PSD Efficiencies vs Cut current"), current_curve)
  end

  curves = EfficiencyCollection[]
  if current_curve != nothing
    push!(curves, current_curve)
  end
  append!(curves, peak_effs)
  push!(curves, random_guessing())

  plot_efficiencies(joinpath(dir, "PSD Efficiency comparison"), curves...)
end

function Plots.plot(curve::EfficiencyCurve)
    plot((curve.cut_values, curve.efficiencies), label=curve.name)
    xaxis!("Cut value")
    yaxis!("Efficiency")
end

function Plots.plot!(curve::EfficiencyCurve)
    plot!(curve.cut_values, curve.efficiencies, label=curve.name)
end


text:
using Test
using AlgebraicDependencies

testnum = []
testdep = []

macro num(numbers...)
    push!(testnum, collect(numbers))
end

macro dep(deps...)
    push!(testdep, collect(deps))
end

gsym(n::Int) = [Symbol("v$i") for i in 1:n]

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

@num 1/2 2
@dep v1*v2-1

@num 1/2 1 2
@dep v2-1 v1*v3-1

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

for (n, d) in zip(testnum, testdep)
    @test dependencies(n; variables=gsym(length(n))) == d
end

ext: jl


text:
# load_constraints.jl
# Gets city information from manhattan.jld object and travel time information from Travel time CSV
# Authored by Arthur J Delarue on 6/8/15

function saveManhattan(pb::TaxiProblem, name::AbstractString; compress=false)
  save("Cities/Saved/$name.jld", "pb", pb, compress=compress)
end

function loadManhattan(name::AbstractString)
  pb = load("Cities/Saved/$name.jld","pb")
  return pb
end

function saveRoadTimes(times::SparseMatrixCSC{Float64,Int}, name::AbstractString; compress=false)
  save("Outputs/$name.jld", "times", times, compress=compress)
end

function loadRoadTimes(name::AbstractString)
  times = load("Outputs/$name.jld","times")
  return times
end

MANHATTAN_NODES = 6134
# Speed limits in km/h from OpenStreetMap/src/speeds.jl
URBAN_SPEEDS = [95,72,48,32,22,12,8,5]

function loadCityGraph(;fromScratch::Bool=false, useShortestPaths::Bool=true,fixSpeeds=false)
	"""
	Load Manhattan object and return it
	"""
	println("**** Loading city graph ****")
	if fromScratch
		manhattan = Manhattan(sp=useShortestPaths)
		# if useShortestPaths
		# 	sp = manhattan.sp
		# else
		# 	sp = parallelShortestPathsAuto(manhattan.network, manhattan.roadTime, manhattan.roadCost)
		# 	manhattan.sp = sp
		# end
		if fixSpeeds
			# Make sure the speeds are in URBAN_SPEEDS
			# If not match speed to closest actual speed
			for i=1:MANHATTAN_NODES,j=1:MANHATTAN_NODES
				if manhattan.distances[i,j] != 0
					# Don't forget to convert from m/s to km/h
					if manhattan.roadTime[i,j] != 0
						speed = 3.6 * manhattan.distances[i,j]/manhattan.roadTime[i,j]
					else
						speed = 3.6 * manhattan.distances[i,j]/0.1
					end
					road_type = indmin(abs(URBAN_SPEEDS - speed))
					if manhattan.roadTime[i,j] == 0
						println(manhattan.roadTime[i,j], "\t", manhattan.distances[i,j], "\t", speed, "\t", URBAN_SPEEDS[road_type])
					end
					manhattan.roadTime[i,j] = manhattan.distances[i,j]/URBAN_SPEEDS[road_type]
				end
			end
			shortestPaths!(manhattan)
		end
		saveManhattan(manhattan, "full_manhattan")
	else
		manhattan = loadManhattan("full_manhattan")
	end
	return manhattan
end

function loadTravelTimeData(;radius::Int=40, times::AbstractString= "1214", min_rides::Int=4, num_nodes::Int=MANHATTAN_NODES, preprocess::Bool=false, num_clusters::Int64 = 50, sampleSize::Int64=10000)
	# Hardcoded for simplicity. Map of Manhattan not likely to change anytime soon
	"""
	Return array of travel times between pairs of nodes A,B obtained from NYC dataset
	"""
	println("**** Loading travel times ****")
	# Load Manhattan highway nodes
	highwayNodes = load("Cities/Saved/highwayNodes.jld", "highwayNodes")
	if preprocess
		dataFile = "../Travel_times/travel_times_r$(radius)_wd_$(times)_clust$(num_clusters)_rides$(sampleSize)_minr$(min_rides).csv"
	else
		dataFile = "../Travel_times/travel_times_r$(radius)_wd_$(times).csv"
	end
	println("-- Loading from $dataFile...")
	data = readcsv(dataFile)
	nodePairs = data[:,1]
	average = data[:,2]
	minRides = data[:,5]

	# Put travel time data in an array
	travelTimes = zeros(num_nodes,num_nodes)
	for (i, element) in enumerate(nodePairs)
		node = split(element, ";")
		# Remember +1 because Python indexes from 0 but Julia from 1
		src = int(node[1]) + 1
		dest = int(node[2]) + 1
		if minRides[i] >= min_rides && average[i] >= 60 && src != dest && !(src in highwayNodes) && !(dest in highwayNodes)
			travelTimes[src, dest] = average[i]
		end
	end
	return travelTimes
end

function loadNewTravelTimeData(;trainOrTest::AbstractString="training", radius::Int=40, times::AbstractString="1214", num_nodes::Int=MANHATTAN_NODES, preprocess::Bool=false, num_clusters::Int64=50, sampleSize::Int64=10000, min_rides::Int=4, saveTestingMatrix::Bool=false, loadTestingMatrixDirectly=false)
	"""
	Return array of travel times between pairs of nodes A,B obtained from NYC dataset, as well as number of rides for each node pair.
	"""
	println("**** Loading $(trainOrTest) data ****")
	# Load Manhattan highway nodes
	highwayNodes = load("Cities/Saved/highwayNodes.jld", "highwayNodes")
	if preprocess
		dataFile = "../Travel_times/$(trainOrTest)_r$(radius)_wd_$(times)_clust$(num_clusters)_rides$(sampleSize)_minr$(min_rides).csv"
	else
		dataFile = "../Travel_times/$(trainOrTest)_r$(radius)_wd_$(times).csv"
	end
	# Special case: load testing matrix directly
	if preprocess == false && loadTestingMatrixDirectly == true
		println("-- Loading directly from JLD archive")
		travelTimes = load("../Travel_times/$(trainOrTest)_r$(radius)_wd_$(times).jld", "travelTimes")
		numRides = load("../Travel_times/$(trainOrTest)_r$(radius)_wd_$(times).jld", "numRides")
		return travelTimes, numRides
	end
	println("-- Loading from $dataFile...")
	data = readcsv(dataFile)
	nodePairs = data[:,1]
	averageTime = data[:,2]
	minRides = data[:,5]
	travelTimes = zeros(num_nodes,num_nodes)
	numRides = zeros(Int, (num_nodes, num_nodes))
	for (i, element) in enumerate(nodePairs)
		node = split(element, ";")
		src = int(node[1]) + 1
		dest = int(node[2]) + 1
		if src != dest && !(src in highwayNodes) && !(dest in highwayNodes) && averageTime[i] >= 60
			travelTimes[src, dest] = averageTime[i]
			numRides[src, dest] = minRides[i]
		end
	end
	if preprocess == false && saveTestingMatrix == true
		println("-- Saving to JLD file for speed")
		save("../Travel_times/$(trainOrTest)_r$(radius)_wd_$(times).jld", "travelTimes", travelTimes, "numRides", numRides)
	end
	return travelTimes, numRides
end

function chooseConstraints(travelTimes::Array{Float64, 2}, numRides::Array{Int, 2}; num_nodes::Int = MANHATTAN_NODES, sample_size::Int = 5000)
	"""
	Given a matrix of travel times and number of rides, randomly selects sample_size nonzero entries
	"""
	newTravelTimes = zeros(num_nodes, num_nodes)
	newNumRides = zeros(Int, (num_nodes, num_nodes))
	for i = 1:sample_size
		startNode = rand(1:num_nodes)
		endNode = rand(1:num_nodes)
		while true
			# Keep them if not added yet
			if numRides[startNode, endNode] > 0 && newNumRides[startNode, endNode] == 0
				break
			end
			startNode = rand(1:num_nodes)
			endNode = rand(1:num_nodes)
		end
		newNumRides[startNode, endNode] = numRides[startNode, endNode]
		newTravelTimes[startNode, endNode] = travelTimes[startNode, endNode]
	end
	return newTravelTimes, newNumRides
end

max_stars_repo_name: sebmart/NYC-timing

max_issues_repo_name: sebmart/NYC-timing

max_forks_repo_name: sebmart/NYC-timing


text:
## Julia program for (Envelope) Accept-Reject Method
## author: weiya <szcfweiya@gmail.com>
## date: 2018-08-19

function AccRej(f::Function, g::Function, M)
    ## use normal distribution N(0, 1) as sampling function g
    while true
        x = randn()
        u = rand()
        cutpoint = f(x)/(M*g(x))
        if u <= cutpoint
            return(x)
        end
    end
end

function EnvAccRej(f::function, M, gl::function)
    while true
        x = randn() # still assume gm is N(0,1)
        u = rand()
        cutpoint1 = gl(x)/(M*g(x))
        cutpoint2 = f(x)/(M*g(x))
        if u <= cutpoint1
            return(x)
        elseif u <= cutpoint2
            return(x)
        end
    end
end

## density function of N(0, 1)
function g(x)
    return(exp(-0.5*x^2)/sqrt(2*pi))
end

## example function and ignore the normalized constant
function f(x)
    return(exp(-x^2/2)*(sin(6*x)^2 + 3*cos(x)^2*sin(4*x)^2 + 1))
end

## example
N = 500;
data = ones(N);
for i = 1:500
    data[i] = AccRej(f, g, sqrt(2*pi)*5)
end

max_stars_repo_name: szcf-weiya/MCMC

max_issues_repo_name: szcf-weiya/MCMC

max_forks_repo_name: szcf-weiya/MCMC


text:
include("./deeper-cut/deeper-cut.config.jl")

function trainresults(
    save_tag,
    model,
    data_trn,
    data_tst,
    learning_rate_per_epoch,
    optimizer = momentum,
    should_save = false;
    accuracy_func = accuracy,
    error_func,
    secondary_accuracy_func,
    data_items_trn,
    data_items_val,
    check_loss = true
)
    GC.gc(true)
    epoch_count = 0
    # TODO: we can make better use of it no need to run forward pass multiple times
    function snapshot(i)

        trn_loss = check_loss ? model(data_trn) : 0
        tst_loss = check_loss ? model(data_trn) : 0

        naive_trn_accuracy = compute_accuracy_in_training(model, data_trn, accuracy_func)
        naive_tst_accuracy = compute_accuracy_in_training(model, data_tst, accuracy_func)

        pck_trn_accuracy, pck_trn_distribution =
            data_items_trn == nothing ? (0, 0) :
            compute_secondary_accuracy_in_training(
                model,
                data_trn,
                data_items_trn,
                secondary_accuracy_func,
            )

        pck_val_accuracy, pck_val_distribution = data_items_val == nothing ?  (0, 0) : compute_secondary_accuracy_in_training(
            model,
            data_tst,
            data_items_val,
            secondary_accuracy_func,
        )

        if pck_val_distribution != 0
                pck_val_distribution = read_accuracy_results(pck_val_distribution)
        end

        open("$(pre_full_path)$(save_tag)-training_snapshots.txt", "a") do io
            write(io, "***** epoch: $(epoch_count) ***** \n")
            write(io, "trn_loss: $(trn_loss) \n")
            write(io, "tst_loss: $(tst_loss) \n")
            write(io, "naive_trn_accuracy: $(naive_trn_accuracy) \n")
            write(io, "naive_tst_accuracy: $(naive_tst_accuracy)  \n")
            write(io, "pck_val_accuracy: $(pck_val_accuracy)  \n")
            write(io, "$(pck_val_distribution)")
            write(io, "\n")
        end

        epoch_count += 1

        Knet.save("$(pre_full_path)$(save_tag)-training_model.jld2", "$(save_tag)-training_model", model)

        snap_res = (
            0,
            trn_loss,
            tst_loss,
            0,
            0,
            naive_trn_accuracy,
            naive_tst_accuracy,
            pck_trn_accuracy,
            pck_val_accuracy,
        )
        println(snap_res)
        return snap_res
    end
    # compute_error_in_training(model, data_trn, error_func),
    # compute_error_in_training(model, data_tst, error_func),

    results = []
    for (lr, for_epoch) in learning_rate_per_epoch

        training = optimizer(model, ncycle(data_trn, for_epoch), lr = lr,
            gamma = momentum_gamma)

        #   (snapshot(x) for x in takenth(progress(training), length(data_trn))) |> collect

        snapshots = (snapshot(x) for x in takenth(progress(training), length(data_trn)))
        intermediate_res = reshape(collect(Iterators.flatten(snapshots)), (9, :))
        if isempty(results)
            results = intermediate_res
        else
            results = hcat(results, intermediate_res)
        end

    end

    if (should_save)
        Knet.save("$(save_tag)_results.jld2", "$(save_tag)_results", results)
    end
    return results

end

# TODO: modify this for our purposes
function alternative_train!(model, trn, dev, tst...; iteration_count = 1000)
    bestmodel, bestloss = deepcopy(model), loss(model, dev)
    progress!(adam(model, trn), steps = iteration_count) do y
        losses = [loss(model, d) for d in (dev, tst...)]
        if losses[1] < bestloss
            bestmodel, bestloss = deepcopy(model), losses[1]
        end
        return (losses...,)
    end
    return bestmodel
end

function compute_accuracy_in_training(model, data, accuracy_func)
    if accuracy_func !== nothing
        return accuracy_func(model, data)
    else
        return 0
        train_accuracy = 0
    end
end

function compute_secondary_accuracy_in_training(model, data, data_items, accuracy_func)
    if accuracy_func !== nothing
        return accuracy_func(model, data, data_items)
    else
        return (0, 0)
    end
end

function compute_error_in_training(model, data, error_func)
    if error_func !== nothing
        return error_func(model, data)
    else
        return 0
    end
end


text:
using ChaosTools, Test

println("\nTesting expansion entropy...")
@testset "Expansion entropy" begin

@testset "maximalexpansion" begin
@test ChaosTools.maximalexpansion(
    [1. 0. 0. 0. 2.; 0. 0. 3. 0. 0.; 0. 0. 0. 0. 0.; 0. 2. 0. 0. 0.]
    ) ≈ 3.0 * 2.23606797749979 * 2.0
@test ChaosTools.maximalexpansion([1 0; 0 1]) == 1
end

@testset "discrete 1D" begin
    tent_eom(x, p, n) = (x < -0.2 ? -0.6 : (x < 0.4 ? 3x : 2(1-x)))
    tent_jacob(x, p, n) = (x < -0.2 ? 0.0 : (x < 0.4 ? 3.0 : -2.0))
    tent = DiscreteDynamicalSystem(tent_eom, 0.2, nothing, tent_jacob)
    t, tent_meanlist, tent_stdlist = expansionentropy_batch(tent, rand, x -> 0 < x < 1; batches=100, N=100000, steps=60)

    exact_ee = log(2)
    for (i, mean) in enumerate(tent_meanlist)
        @test abs(mean/i - exact_ee) < 0.05
    end

    ee = expansionentropy(tent, rand, x->0<x<1; batches=1000, N=10000, steps=60)
    @test abs(ee - exact_ee) ≤ 0.1
end

@testset "discrete 1D nonchaotic" begin
    expand2_eom(x, p, n) = 2*x
    expand2_jacob(x, p, n) = 2
    expand2 = DiscreteDynamicalSystem(expand2_eom, 0.2, nothing, expand2_jacob)
    _, expand_meanlist, _ = expansionentropy_batch(expand2, rand, x -> 0 < x < 1; batches=100, N=100000, steps=10)

    for (i, mean) in enumerate(expand_meanlist)
        @test -0.1< mean/i < 0.1
    end
end

@testset "discrete 2D" begin
    cat = DynamicalSystemsBase.Systems.arnoldcat()
    cat_gen() = [rand(), rand()]
    cat_inside(x) = true

    _, cat_meanlist, _ = expansionentropy_batch(cat, cat_gen, cat_inside; batches=100, N=100, steps=30, Δt=1)

    exact_ee = log( 1/2 * (3 + sqrt(5)))
    for i in 1:length(cat_meanlist)
        @test abs(cat_meanlist[i]/i - exact_ee) < 1e-6
    end

    ee = expansionentropy(cat, cat_gen, cat_inside; batches=100, N=100, steps=30, Δt=1)
    @test ee ≈ exact_ee
end

@testset "continuous chaotic (lorenz)" begin
    lor = DynamicalSystemsBase.Systems.lorenz()
    lor_gen() = [rand()*40-20, rand()*60-30, rand()*50]
    lor_isinside(x) = -20 < x[1] < 20 && -30 < x[2] < 30 && 0 < x[3] < 50
    ee = expansionentropy(lor, lor_gen, lor_isinside; batches=100, N=100, steps=40, Δt=1.0)

    @test abs(ee - 0.9) < 0.07
end

@testset "continuous regular (henon)" begin
    hh = Systems.henonheiles([0.0, 0.1, 0.5, 0.0])
    tr = trajectory(hh, 100.0, Δt = 0.1, Ttr=200)
    x, y, z, w = columns(tr)
    gen, restr = boxregion(map(minimum, columns(tr)), map(maximum, columns(tr)))

    t, m, s = expansionentropy_batch(hh, gen, restr;
              batches=1000, N=100, steps=40, Δt=0.1, Ttr=1.0)
    @test all(m[10:20] .< 0.01)
end
end


text:
"""
TableSchema module
https://github.com/frictionlessdata/tableschema-jl
"""
module TableSchema

export Table
export Schema
export Field

export TableValidationException
export ConstraintError
export SchemaError

using Base.Iterators: filter
using Base.Iterators: Repeated, repeated

using JSON
using Dates

import HTTP: request
import DelimitedFiles: readdlm, writedlm

const DEFAULT_TYPE = "string"
const DEFAULT_FORMAT = "default"
const MAX_ROWS_INFER = 100

include("constraints.jl")
include("field.jl")
include("exceptions.jl")
include("schema.jl")
include("table.jl")
include("validators.jl")

end # module


text:
# todo = parse(Int64, ARGS[1]):10:290
todo = 280:280

for doing ∈ todo
    println(doing)
    try
        run(`julia aespa.jl $doing`)
    catch LoadError
        open("../../바탕 화면/error $doing.txt", "w")
    end
end

max_stars_repo_name: rmsmsgood/aespa

max_issues_repo_name: rmsmsgood/aespa

max_forks_repo_name: rmsmsgood/aespa


text:
module DistributionsAD

using PDMats,
      ForwardDiff,
      LinearAlgebra,
      Distributions,
      Random,
      Combinatorics,
      SpecialFunctions,
      StatsFuns,
      Compat,
      Requires

using Tracker: Tracker, TrackedReal, TrackedVector, TrackedMatrix, TrackedArray,
                TrackedVecOrMat, track, @grad, data
using SpecialFunctions: logabsgamma, digamma
using LinearAlgebra: copytri!, AbstractTriangular
using Distributions: AbstractMvLogNormal,
                     ContinuousMultivariateDistribution
using DiffRules, SpecialFunctions, FillArrays
using ForwardDiff: @define_binary_dual_op # Needed for `eval`ing diffrules here
using Base.Iterators: drop

import StatsFuns: logsumexp,
                  binomlogpdf,
                  nbinomlogpdf,
                  poislogpdf,
                  nbetalogpdf
import Distributions: MvNormal,
                      MvLogNormal,
                      poissonbinomial_pdf_fft,
                      logpdf,
                      quantile,
                      PoissonBinomial,
                      Binomial,
                      BetaBinomial,
                      Erlang
import ZygoteRules

export TuringScalMvNormal,
       TuringDiagMvNormal,
       TuringDenseMvNormal,
       TuringMvLogNormal,
       TuringPoissonBinomial,
       TuringWishart,
       TuringInverseWishart,
       arraydist,
       filldist

include("common.jl")
include("univariate.jl")
include("multivariate.jl")
include("mvcategorical.jl")
include("matrixvariate.jl")
include("flatten.jl")
include("arraydist.jl")
include("filldist.jl")
@init @require ReverseDiff = "37e2e3b7-166d-5795-8a7a-e32c996b4267" begin
    include("reversediff.jl")
end

end


text:
# Print a floating point number at fixed precision. Pretty much equivalent to
# @sprintf("%0.$(precision)f", x), without the macro issues.
function format_fixed(x::FloatingPoint, precision::Integer)
    @assert precision >= 0

    if x == Inf
        return "∞"
    elseif x == -Inf
        return "-∞"
    elseif isnan(x)
        return "NaN"
    end

    Base.Grisu.@grisu_ccall x Base.Grisu.FIXED precision
    point, len, digits = (Base.Grisu.POINT[1], Base.Grisu.LEN[1], Base.Grisu.DIGITS)
    ss = Uint8[]
    if x < 0
        push!(ss, '-')
    end

    append!(ss, digits[1:min(point, len)])

    if point > len
        for _ in len:point-1
            push!(ss, '0')
        end
    elseif point < len
        if point <= 0
            push!(ss, '0')
        end
        push!(ss, '.')
        if point < 0
            for _ in 1:-point
                push!(ss, '0')
            end
            append!(ss, digits[1:len])
        else
            append!(ss, digits[point+1:len])
        end
    end

    trailing_zeros = precision - max(0, len - point)
    if trailing_zeros > 0 && point >= len
        push!(ss, '.')
    end

    for _ in 1:trailing_zeros
        push!(ss, '0')
    end

    bytestring(ss)
end

# Print a floating point number in scientific notation at fixed precision. Sort of equivalent
# to @sprintf("%0.$(precision)e", x), but prettier printing.
function format_fixed_scientific(x::FloatingPoint, precision::Integer)
    if x == 0.0
        return "0"
    elseif x == Inf
        return "∞"
    elseif x == -Inf
        return "-∞"
    elseif isnan(x)
        return "NaN"
    end

    mag = log10(abs(x))
    if mag < 0
        grisu_precision = precision + abs(iround(mag))
    else
        grisu_precision = precision
    end

    Base.Grisu.@grisu_ccall x Base.Grisu.FIXED grisu_precision
    point, len, digits = (Base.Grisu.POINT[1], Base.Grisu.LEN[1], Base.Grisu.DIGITS)

    @assert len > 0

    ss = Char[]
    if x < 0
        push!(ss, '-')
    end

    push!(ss, digits[1])
    if precision > 1
        push!(ss, '.')
    end

    for i in 2:len
        push!(ss, digits[i])
    end

    for i in (len+1):precision
        push!(ss, '0')
    end

    string(utf8(ss), "×10<sup>$(point - 1)</sup>")
end

formatter() = string

# Nicely format an print some numbers.
function formatter{T<:FloatingPoint}(xs::AbstractArray{T}; fmt=:auto)
    if length(xs) <= 1
        return string
    end

    # figure out the lowest suitable precision
    delta = Inf
    finite_xs = filter(isfinite, xs)
    for (x0, x1) in zip(finite_xs, drop(finite_xs, 1))
        delta = min(x1 - x0, delta)
    end
    x_min, x_max = concrete_minimum(xs), concrete_maximum(xs)

    x_min, x_max, delta = (float64(float32(x_min)), float64(float32(x_max)),
        float64(float32(delta)))

    if !isfinite(x_min) || !isfinite(x_max) || !isfinite(delta)
        error("At least one finite value must be provided to formatter.")
    end

    if fmt == :auto
        if abs(log10(x_max - x_min)) > 4
            fmt = :scientific
        else
            fmt = :plain
        end
    end

    if fmt == :plain
        # SHORTEST_SINGLE rather than SHORTEST to crudely round away tiny innacuracies
        Base.Grisu.@grisu_ccall delta Base.Grisu.SHORTEST_SINGLE 0
        precision = max(0, Base.Grisu.LEN[1] - Base.Grisu.POINT[1])

        return x -> format_fixed(x, precision)
    elseif fmt == :scientific
        Base.Grisu.@grisu_ccall delta Base.Grisu.SHORTEST_SINGLE 0
        delta_magnitude = Base.Grisu.POINT[1]

        Base.Grisu.@grisu_ccall x_max Base.Grisu.SHORTEST_SINGLE 0
        x_max_magnitude = Base.Grisu.POINT[1]

        precision = 1 + max(0, x_max_magnitude - delta_magnitude)

        return x -> format_fixed_scientific(x, precision)
    end
end

# Print dates
function formatter{T<:Date}(ds::AbstractArray{T}; fmt=nothing)
    if isa(fmt, Function)
        return fmt
    end

    const month_names = [
        "January", "February", "March", "April", "May", "June", "July",
        "August", "September", "October", "November", "December"
    ]

    const month_abbrevs = [
        "Jan", "Feb", "Mar", "Apr", "May", "June",
        "July", "Aug", "Sept", "Oct", "Nev", "Dec"
    ]

    day_all_1   = all(map(d -> day(d) == 1, ds))
    month_all_1 = all(map(d -> month(d) == 1, ds))
    years = Set()
    for d in ds
        push!(years, year(d))
    end

    if day_all_1 && month_all_1
        # only label years
        function format(d)
            buf = IOBuffer()
            print(buf, year(d))
            takebuf_string(buf)
        end
    elseif day_all_1
        # label months and years
        function format(d)
            buf = IOBuffer()
            if d == ds[1] || month(d) == 1
                print(buf, month_abbrevs[month(d)], " ", year(d))
            else
                print(buf, month_abbrevs[month(d)])
            end
            takebuf_string(buf)
        end
    else
        function format(d)
            buf = IOBuffer()
            if d == ds[1] || (month(d) == 1 && day(d) == 1)
                print(buf, month_abbrevs[month(d)], " ", day(d), " ", year(d))
            elseif day(d) == 1
                print(buf, month_abbrevs[month(d)], " ", day(d))
            else
                print(buf, day(d))
            end
            takebuf_string(buf)
        end
    end
end

# Catchall
function formatter(xs::AbstractArray; fmt=nothing)
    function format(x)
        buf = IOBuffer()
        print(buf, x)
        takebuf_string(buf)
    end

    format
end


text:
# ---
# title: 14. Longest Common Prefix
# id: problem14
# author: Tian Jun
# date: 2020-10-31
# difficulty: Easy
# categories: String
# link: <https://leetcode.com/problems/longest-common-prefix/description/>
# hidden: true
# ---
#
# Write a function to find the longest common prefix string amongst an array of
# strings.
#
# If there is no common prefix, return an empty string `""`.
#
#
#
# **Example 1:**
#
#
#
#     Input: strs = ["flower","flow","flight"]
#     Output: "fl"
#
#
# **Example 2:**
#
#
#
#     Input: strs = ["dog","racecar","car"]
#     Output: ""
#     Explanation: There is no common prefix among the input strings.
#
#
#
#
# **Constraints:**
#
#   * `0 <= strs.length <= 200`
#   * `0 <= strs[i].length <= 200`
#   * `strs[i]` consists of only lower-case English letters.
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end


text:
using BinDeps

@BinDeps.setup

libspg = library_dependency("libsymspg")
provides(Sources, URI("https://sourceforge.net/projects/spglib/files/spglib/spglib-1.9/spglib-1.9.4.tar.gz"), libspg)
provides(BuildProcess, Autotools(libtarget = joinpath("src", "libsymspg.la")), libspg)
@BinDeps.install Dict(:libsymspg => :spglib)


text:
function write(obj::EphemerisFileWriter, arg0::Appendable, arg1::EphemerisFile)
    return jcall(obj, "write", void, (Appendable, EphemerisFile), arg0, arg1)
end

function write(obj::EphemerisFileWriter, arg0::JString, arg1::EphemerisFile)
    return jcall(obj, "write", void, (JString, EphemerisFile), arg0, arg1)
end


text:
function solve(filename)
    values = readlines(filename*".txt")
    values = split.(values, " | ")

    input = String[]
    output = String[]
    for line in values
        input = push!(input, line[1])
        output = push!(output, line[2])
    end
    input = split.(input, ' ')
    output = split.(output, ' ')

    count = 0
    for line in output
        for el in line
            if length(el) == 2 || length(el) == 4 || length(el) == 3 || length(el) == 7
                count += 1
            end
        end
    end
    println("Count of 1, 4, 7, 8: "*string(count))
end

function find_patterns(input)
    codesperline = []
    zero = one = two = three = four = five = six = seven = eight = nine = ""
    upper_line = mid_line = lower_line = right_upper = right_lower = left_upper = left_lower = ""
    for line in input
        #one
        for el in line
            if length(el) == 2
                one = join(sort(collect(el)))
                break
            end
        end
        #seven
        for el in line
            if length(el) == 3
                seven = join(sort(collect(el)))
                upper_line = replace(el, one => "")
                break
            end
        end
        #four
        for el in line
            if length(el) == 4
                four = join(sort(collect(el)))
                break
            end
        end
        #eight
        for el in line
            if length(el) == 7
                eight = join(sort(collect(el)))
                break
            end
        end
        #three is length five and contains one
        length_five_candidates = []
        length_six_candidates = []
        for el in line
            if length(el) == 5
                length_five_candidates = push!(length_five_candidates, join(sort(collect(el))))
            elseif length(el) == 6
                length_six_candidates = push!(length_six_candidates, join(sort(collect(el))))
            end
        end
        #three
        for el in length_five_candidates
            if stringcontains(one, el)
                three = el
                setdiff!(length_five_candidates, three)
                break
            end
        end
        #six
        for el in length_six_candidates
            if stringcontains(one, el) == false
                six = el
                setdiff!(length_six_candidates, six)
                break
            end
        end
        #nine
        for el in length_six_candidates
            if stringcontains(three, el)
                nine = el
                setdiff!(length_six_candidates, nine)
                break
            end
        end
        #zero
        for el in length_six_candidates
            if stringcontains(nine, el) == false && stringcontains(six, el) == false
                zero = el
                break
            end
        end

        #two
        for char in nine
            if occursin(char, three) == false
                left_upper = char
                break
            end
        end

        for el in length_five_candidates
            if stringcontains(left_upper, el) == false && el != three
                two = el
                setdiff!(length_five_candidates, two)
                break
            end
        end
        #five
        for el in length_five_candidates
            if stringcontains(three, el) == false && stringcontains(two, el) == false
                five = el
                break
            end
        end

        numbers = [zero, one, two, three, four, five, six, seven, eight, nine]
        for (i,n) in enumerate(numbers)
            numbers[i] = join(sort(collect(n)))
        end
        codesperline = push!(codesperline, numbers)
    end
    return codesperline
end

function stringcontains(ptrn, str)
    res = true
    for char in ptrn
        if occursin(char, str) == false
            res = false
            break
        end
    end
    return res
end

function decode(output, codesperline)
    output_numbers = []
    for (i,line) in enumerate(output)
        numbers = codesperline[i]
        decoded = ""
        for el in line
            el = join(sort(collect(el)))
            for i in 0:9
                if el == numbers[i+1]
                    decoded = decoded * string(i)
                    break
                end
            end
        end
        decoded_int = parse(Int, decoded)
        output_numbers = push!(output_numbers, decoded_int)
    end

    sum = 0
    for number in output_numbers
        sum += number
    end
    return sum
end

function solve2(filename)
    values = readlines(filename*".txt")
    values = split.(values, " | ")

    input = String[]
    output = String[]
    for line in values
        input = push!(input, line[1])
        output = push!(output, line[2])
    end
    input = split.(input, ' ')
    output = split.(output, ' ')

    codes = find_patterns(input)
    res = decode(output, codes)
    println("Sum of all numbers: "*string(res))
end

solve("smalltest")
solve("input")

solve2("smalltest")
solve2("input")


text:
struct NoisyOracle <: Oracle
    oracle::Oracle
    p::Float64
end

function NoisyOracle(data, labels, params::Dict{Symbol, Any})
    suboracle = initialize_oracle(eval(params[:subtype]), data, labels, get(params, :subparams, Dict{Symbol, Any}()))
    return NoisyOracle(suboracle, params[:p])
end

function ask_oracle(oracle::NoisyOracle, queries::Union{Vector{Int}, Array{T, 2}})::Vector{Symbol} where T <: Real
    labels = ask_oracle(oracle.oracle, queries)
    flipped_labels = [x == :inlier ? :outlier : :inlier for x in labels]
    mask = rand(length(labels)) .< oracle.p
    labels[mask] .= flipped_labels[mask]
    return labels
end


text:
@testset "Pruning functions tests" begin
    Random.seed!(0)
    X = DualOptimalFiltering.generate_CIR_trajectory(range(0, stop = 2, length = 20), 3, 3., 0.5, 1);
    Y = map(λ -> rand(Poisson(λ),10), X);
    data = Dict(zip(range(0, stop = 2, length = 20), Y))
    Λ_of_t, wms_of_t, θ_of_t = DualOptimalFiltering.filter_CIR(3., 0.5, 1.,1.,data);

    times = Λ_of_t |> keys |> collect |> sort
    Λ_of_t_kept, wms_of_t_kept = DualOptimalFiltering.keep_above_threshold(Λ_of_t[times[end]], wms_of_t[times[end]], 0.001)
    @test minimum(wms_of_t_kept) >= 0.001
    Λ_of_t_kept, wms_of_t_kept = DualOptimalFiltering.keep_fixed_number_of_weights(Λ_of_t[times[end]], wms_of_t[times[end]], 200)
    @test length(wms_of_t_kept |> unique) <= 200
    Λ_of_t_kept, wms_of_t_kept = DualOptimalFiltering.keep_fixed_number_of_weights(Λ_of_t[times[end]], wms_of_t[times[end]], 2)
    @test length(wms_of_t_kept |> unique) <= 2
    Λ_of_t_kept, wms_of_t_kept = DualOptimalFiltering.keep_fixed_number_of_weights(Λ_of_t[times[end]], wms_of_t[times[end]], 1)
    @test length(wms_of_t_kept |> unique) <= 1
    Λ_of_t_kept, wms_of_t_kept = DualOptimalFiltering.keep_fixed_fraction(Λ_of_t[times[end]], wms_of_t[times[end]], 0.95)
    @test sum(wms_of_t_kept) >= 0.95
    Λ_of_t_kept, wms_of_t_kept = DualOptimalFiltering.keep_fixed_fraction(Λ_of_t[times[end]], wms_of_t[times[end]], 0.99)
    @test sum(wms_of_t_kept) >= 0.99

end;

max_stars_repo_name: konkam/DualOptimalFiltering

max_issues_repo_name: konkam/DualOptimalFiltering

max_forks_repo_name: konkam/DualOptimalFiltering


text:
#__precompile__()
module SDEs

include("noises.jl")
include("callbacks.jl")
include("solvers.jl")
include("explicitEM.jl")
include("pcEuler.jl")

export Solver
export EM
export PCE
export default_seed
export integrate

end # module

max_stars_repo_name: onoderat/SDEs

max_issues_repo_name: onoderat/SDEs

max_forks_repo_name: onoderat/SDEs


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

using Test

struct MyMaterial <: AbstractMaterial
    E :: Float64
end

function MyMaterial()
    error("MyMaterial needs at least `E` defined.")
end

mat = Material(MyMaterial, (E=210.0,))
@test mat.properties.E == 210.0


text:
cfg = NeuroEvolution.get_config("../cfg/hyperneat.yaml")

@testset "HyperNEAT Layer" begin
    cfg = NeuroEvolution.get_config("../cfg/hyperneat.yaml")
    l = Layer(10, 0., cfg["hn_activ_func"])
    @test l.size == 10
    @test length(l.values)==10
    @test length(l.biases)==10
    @test l.depth == 0
end

@testset "HyperNEAT Connection" begin
    # Empty layer
    l1 = Layer(10, 0., cfg["hn_activ_func"])
    l2 = Layer(8, 1., cfg["hn_activ_func"])
    c = Connection(l1, l2)
    @test c.layer_from == l1
    @test c.layer_to == l2
    @test size(c.weights) == (10, 8)

    # Generated layer
    f = x -> rand()
    c = Connection(f, l1, l2)
    @test c.layer_from == l1
    @test c.layer_to == l2
    @test size(c.weights) == (10, 8)
end

@testset "HyperNEAT Network" begin
    # All possible connections
    cfg["hn_link_all_layers"]= true
    net = GridNetwork(cfg)
    @test length(net.layers)==length(cfg["hn_hidden_layers"])+2
    @test length(net.connections)==factorial(length(cfg["hn_hidden_layers"])+1)
    X = zeros(cfg["n_in"]) .+ 1
    y = process(net, X)
    @test length(y) == cfg["n_out"]
    @test !all(y .== 0.)
    for c in net.connections
        @test c.layer_from.depth < c.layer_to.depth
    end

    # Only feed-forward
    cfg["hn_link_all_layers"]= false
    net = GridNetwork(cfg)
    @test length(net.layers)==length(cfg["hn_hidden_layers"])+2
    @test length(net.connections)==length(cfg["hn_hidden_layers"])+1
    X = zeros(cfg["n_in"]) .+ 1
    y = process(net, X)
    @test length(y) == cfg["n_out"]
    @test !all(y .== 0.)
    for c in net.connections
        @test c.layer_from.depth < c.layer_to.depth
    end
end

@testset "HyperNEAT Evaluate" begin
    e::Evolution = HyperNEAT(cfg, fitness_xor, 2)
    e.evaluate(e)
    @test all(getfield.(e.population, :fitness) .!= [-Inf])
end

@testset "HyperNEAT run" begin
    cfg = NeuroEvolution.get_config("../cfg/hyperneat.yaml")

    e = HyperNEAT(cfg, fitness_xor, 2)
    Cambrian.run!(e)

    X, y = xor_dataset(2, 100)
    max_f = sum(log_fitness.(y, y)) / length(X)
    y_false = []
    for i in y
        push!(y_false, 1 .- i)
    end
    min_f = sum(log_fitness.(y, y_false)) / length(X)
    if cfg["verbose"]
        println("\n-- HyperNEAT on XOR --")
        println("Worst possible fitness: ", min_f)
        println("Best  possible fitness: ", max_f)
    end
    # Analyse results
    best = sort(e.population, rev=true)
    @test best[1].fitness[1] <= max_f
    @test best[1].fitness[1] > min_f

    y = process(best[1], [1., 1.])
    @test typeof(y[1])==Float64

    if cfg["verbose"]
        println("Final fitness: ", best[1].fitness[1])
        println("Genome size: ", length(best[1].genes))
        println("Max innovation: ", e.cfg["innovation_max"])
        println("Species: ", length(e.cfg["Species"]))
        println("Best(0, 0)= ", y[1])
        println("\nBest individuals: ", getfield.(best[1:5], :fitness))
        println("\nConnections in Best: ")
        for g in values(best[1].genes)
            println(g.origin, " -- ", g.weight, " -> ", g.destination)

        end
    end
end

@testset "Cambrian.GA x HyperNEAT" begin
    cfg = NeuroEvolution.get_config("../cfg/ga.yaml")
    cfg["hyperneat"]=true
    cfg["hn_activ_func"]=sigmoid
    e = GA_NEAT(HyperNEATIndividual, cfg, fitness_xor; id="test")
    step!(e)
    @test length(e.population) == cfg["n_population"]
    run!(e)
    best = sort(e.population, rev=true)

    X, y = xor_dataset(cfg["n_in"], 100)
    max_f = sum(log_fitness.(y, y)) / length(X)
    y_false = []
    for i in y
        push!(y_false, 1 .- i)
    end
    min_f = sum(log_fitness.(y, y_false)) / length(X)
    @test best[1].fitness[1] <= max_f
    @test best[1].fitness[1] > min_f

    y = process(best[1], [1., 1.])
    @test typeof(y[1])==Float64
end


text:
module Day12

moons = Dict(enumerate([[-13, 14, -7], [-18, 9, 0], [0, -3, -3], [-15, 3, -13]]))

potenergy(moon) = sum(abs, moon)
kinenergy(vel) = sum(abs, vel)
totalenergy(m, v) = sum(potenergy.(values(m)) .* kinenergy.(values(v)))

function step!(moons, velocities)
    for i in 1:4, j in 1:4
        i == j && continue
        velocities[i] .+= Int.(map(sign, moons[j] .- moons[i]))
    end

    for k in keys(moons)
        moons[k] .+= velocities[k]
    end
end

function solvea(moons)
    velocities = Dict(i => [0, 0, 0] for i in 1:4)
    for _ in 1:1000
        step!(moons, velocities)
    end

    totalenergy(moons, velocities)
end

function solveb(moons)
    velocities = Dict(i => [0, 0, 0] for i in 1:4)
    seen = [Set(), Set(), Set()]
    stop = falses(3)

    while ! all(stop)
        vals = vcat(values(moons)..., values(velocities)...) |> x -> reshape(x, :, 8)
        for axis in 1:3
            if ! stop[axis]
                vals[axis, :] in seen[axis] && (stop[axis] = true)
                push!(seen[axis], vals[axis, :])
            end
        end
        step!(moons, velocities)
    end

    reduce(lcm, length.(seen))
end

end

max_stars_repo_name: AnAverageHuman/competitive

max_issues_repo_name: AnAverageHuman/competitive

max_forks_repo_name: AnAverageHuman/competitive


text:
"""
**Student-T likelihood**

Student-t likelihood for regression: ``\\frac{\\Gamma((\\nu+1)/2)}{\\sqrt{\\nu\\pi}\\Gamma(\\nu/2)}\\left(1+t^2/\\nu\\right)^{(-(\\nu+1)/2)}``
see [wiki page](https://en.wikipedia.org/wiki/Student%27s_t-distribution)

```julia
StudentTLikelihood(ν::T,σ::Real=one(T)) #ν is the number of degrees of freedom
#σ is the variance for local scale of the data.
```

---

For the analytical solution, it is augmented via:
```math
p(y|f,\\omega) = \\mathcal{N}(y|f,\\omega)
```
Where ``\\omega \\sim \\mathcal{IG}(\\frac{\\nu}{2},\\frac{\\nu}{2})`` where ``\\mathcal{IG}`` is the inverse gamma distribution
See paper [Robust Gaussian Process Regression with a Student-t Likelihood](http://www.jmlr.org/papers/volume12/jylanki11a/jylanki11a.pdf)
"""
struct StudentTLikelihood{T<:Real} <: RegressionLikelihood{T}
    ν::T
    α::T
    σ::T
    β::LatentArray{Vector{T}}
    θ::LatentArray{Vector{T}}
    function StudentTLikelihood{T}(ν::T,σ::T=one(T)) where {T<:Real}
        new{T}(ν,(ν+one(T))/2.0,σ)
    end
    function StudentTLikelihood{T}(ν::T,σ::T,β::AbstractVector{<:AbstractVector{T}},θ::AbstractVector{<:AbstractVector{T}}) where {T<:Real}
        new{T}(ν,(ν+one(T))/2.0,σ,β,θ)
    end
end

function StudentTLikelihood(ν::T,σ::T=one(T)) where {T<:Real}
    StudentTLikelihood{T}(ν,σ)
end

function init_likelihood(likelihood::StudentTLikelihood{T},inference::Inference{T},nLatent::Int,nSamplesUsed::Int) where T
    if inference isa AnalyticVI || inference isa GibbsSampling
        StudentTLikelihood{T}(
        likelihood.ν,likelihood.σ,
        [abs2.(T.(rand(T,nSamplesUsed))) for _ in 1:nLatent],
        [zeros(T,nSamplesUsed) for _ in 1:nLatent])
    else
        StudentTLikelihood{T}(likelihood.ν,likelihood.σ)
    end
end

function pdf(l::StudentTLikelihood,y::Real,f::Real)
    tdistpdf(l.ν,(y-f)/l.σ)
end

function Base.show(io::IO,model::StudentTLikelihood{T}) where T
    print(io,"Student-t likelihood")
end

function compute_proba(l::StudentTLikelihood{T},μ::AbstractVector{T},σ²::AbstractVector{T}) where {T<:Real}
    N = length(μ)
    st = TDist(l.ν)
    nSamples = 2000
    μ_pred = zeros(T,N)
    σ²_pred = zeros(T,N)
    temp_array = zeros(T,nSamples)
    for i in 1:N
        # e = expectation(Normal(μ[i],sqrt(σ²[i])))
        # μ_pred[i] = μ[i]
        #
        # σ²_pred[i] = e(x->pdf(LocationScale(x,1.0,st))^2) - e(x->pdf(LocationScale(x,1.0,st)))^2
        if σ²[i] <= 1e-3
            pyf =  LocationScale(μ[i],sqrt(l.σ),st)
            for j in 1:nSamples
                temp_array[j] = rand(pyf)
            end
        else
            d = Normal(μ[i],sqrt(σ²[i]))
            for j in 1:nSamples
                temp_array[j] = rand(LocationScale(rand(d),sqrt(l.σ),st))
            end
        end
        μ_pred[i] = μ[i];
        σ²_pred[i] = cov(temp_array)
    end
    return μ_pred,σ²_pred
end

## Local Updates ##

function local_updates!(model::VGP{<:StudentTLikelihood,<:AnalyticVI})
    model.likelihood.β .= broadcast((Σ,μ,y)->0.5*(Σ+abs2.(μ-y).+model.likelihood.σ*model.likelihood.ν),diag.(model.Σ),model.μ,model.y)
    model.likelihood.θ .= broadcast(β->model.likelihood.α./β,model.likelihood.β)
end

function local_updates!(model::SVGP{<:StudentTLikelihood,<:AnalyticVI})
    model.likelihood.β .= broadcast((K̃,κ,Σ,μ,y)->0.5*(K̃ + opt_diag(κ*Σ,κ) + abs2.(κ*μ-y[model.inference.MBIndices]).+model.likelihood.σ*model.likelihood.ν),model.K̃,model.κ,model.Σ,model.μ,model.y)
    model.likelihood.θ .= broadcast(β->model.likelihood.α./β,model.likelihood.β)
end

function sample_local!(model::VGP{<:StudentTLikelihood,<:GibbsSampling})
    model.likelihood.β .= broadcast((μ::AbstractVector{<:Real},y)->rand.(InverseGamma.(model.likelihood.α,0.5*(abs2.(μ-y).+model.likelihood.σ*model.likelihood.ν))),model.μ,model.y)
    model.likelihood.θ .= broadcast(β->1.0./β,model.likelihood.β)
    return nothing
end

## Global Gradients ##

function cond_mean(model::VGP{<:StudentTLikelihood,<:AnalyticVI},index::Integer)
    return model.likelihood.θ[index].*model.y[index]
end

function ∇μ(model::VGP{<:StudentTLikelihood})
    return hadamard.(model.likelihood.θ,model.y)
end

function cond_mean(model::SVGP{<:StudentTLikelihood,<:AnalyticVI},index::Integer)
    return model.likelihood.θ[index].*model.y[index][model.inference.MBIndices]
end

function ∇μ(model::SVGP{<:StudentTLikelihood,<:AnalyticVI})
    return hadamard.(model.likelihood.θ,getindex.(model.y,[model.inference.MBIndices]))
end

function ∇Σ(model::AbstractGP{<:StudentTLikelihood})
    return model.likelihood.θ
end

function ELBO(model::AbstractGP{<:StudentTLikelihood,<:AnalyticVI})
    return expecLogLikelihood(model) - InverseGammaKL(model) - GaussianKL(model)
end

## ELBO Section ##

function expecLogLikelihood(model::VGP{StudentTLikelihood{T},AnalyticVI{T}}) where T
    tot = -0.5*model.nLatent*model.nSample*log(twoπ)
    tot += -0.5.*sum(broadcast(β->sum(model.nSample*digamma(model.likelihood.α).-log.(β)),model.likelihood.β))
    tot += -0.5.*sum(broadcast((θ,Σ,μ,y)->dot(θ,Σ+abs2.(μ)-2.0*μ.*y-abs2.(y)),model.likelihood.θ,diag.(model.Σ),model.μ,model.y))
    return tot
end

function expecLogLikelihood(model::SVGP{StudentTLikelihood{T},AnalyticVI{T}}) where T
    tot = -0.5*model.nLatent*model.inference.nSamplesUsed*log(twoπ)
    tot += -0.5.*sum(broadcast(β->sum(model.inference.nSamplesUsed*digamma(model.likelihood.α).-log.(β)),model.likelihood.β))
    tot += -0.5.*sum(broadcast((θ,K̃,κ,Σ,κμ,y)->dot(θ,(K̃+opt_diag(κ*Σ,κ)+abs2.(κμ)-2.0*(κμ).*y[model.inference.MBIndices]-abs2.(y[model.inference.MBIndices]))),model.likelihood.θ,model.K̃,model.κ,model.Σ,model.κ.*model.μ,model.y))
    return model.inference.ρ*tot
end

function InverseGammaKL(model::AbstractGP{<:StudentTLikelihood})
    α_p = model.likelihood.ν/2; β_p= α_p*model.likelihood.σ
    model.inference.ρ*sum(broadcast(InverseGammaKL,model.likelihood.α,model.likelihood.β,α_p,β_p))
end

## Numerical Gradients ##

function gradpdf(::StudentTLikelihood,y::Int,f::T) where {T<:Real}
    @error "Not implemented yet"
end

function hessiandiagpdf(::StudentTLikelihood,y::Int,f::T) where {T<:Real}
    @error "Not implemented yet"
end


text:
module QuadTests

using Test
using Jalerkin

qr = TrapezoidQuadrature()
# Access to nodes and weights
@test(get_nodes(qr) == [-1,1])
@test(get_weights(qr) == [1,1])
# Compute $\int_{-1}^1 x dx$ using Trapezoidal quad. rule
@test quad(qr, x->x) ≈ 0  # Exact for P1
# Compute $\int_{-1}^1 x^2 dx$ using Trapezoidal quad. rule
@test !(quad(qr, x->x^2) ≈ 2/3)  # Not exact for P2

# Compute $\int_0^1 x^3 dx$ using Simpson's quad. rule
nodes = [0, 0.5, 1]
weights = [1/6, 2/3, 1/6]
qr = Quadrature{3}(nodes, weights)
@test quad(qr, x->x^3) ≈ 0.25

# Test Gaussian quad. rule on [-1,1]
@test quad(GaussianQuadrature(4), x->x^6) ≈ 2/7

end #module

max_stars_repo_name: rrgalvan/Jalerkin

max_issues_repo_name: rrgalvan/Jalerkin

max_forks_repo_name: rrgalvan/Jalerkin


text:
###
# Cardinalities
###
"""
InfiniteCardinal{K}()

represents the k-th cardinal number.
Note that `InfiniteCardinal <: Integer` to support
being treated as a length in array machinier.
"""
struct InfiniteCardinal{N} <: Integer end

const ℵ₀ = InfiniteCardinal{0}()
const ℵ₁ = InfiniteCardinal{1}()

string(::InfiniteCardinal{0}) = "ℵ₀"
string(::InfiniteCardinal{1}) = "ℵ₁"

show(io::IO, F::InfiniteCardinal{0}) where N = print(io, "ℵ₀")
show(io::IO, F::InfiniteCardinal{1}) where N = print(io, "ℵ₁")

isone(::InfiniteCardinal) = false
iszero(::InfiniteCardinal) = false

signbit(::InfiniteCardinal) = false
sign(::InfiniteCardinal) = 1
angle(::InfiniteCardinal) = 0
abs(a::InfiniteCardinal) = a
zero(::InfiniteCardinal) = 0
one(::Type{<:InfiniteCardinal}) = 1
oneunit(::Type{<:InfiniteCardinal}) = 1
oneunit(::InfiniteCardinal) = 1

isinf(::InfiniteCardinal) = true
isfinite(::InfiniteCardinal) = false

Integer(::Infinity) = InfiniteCardinal{0}()
function Integer(x::RealInfinity)
    signbit(x) && throw(InexactError(:Integer, Integer, x))
    ℵ₀
end
function Integer(x::ComplexInfinity)
    iszero(angle(x)) || throw(InexactError(:Integer, Integer, x))
    ℵ₀
end

==(::InfiniteCardinal{N}, ::InfiniteCardinal{N}) where N = true
==(::InfiniteCardinal, ::InfiniteCardinal)= false
==(::InfiniteCardinal, ::Int) = false
==(::Int, ::InfiniteCardinal) = false
==(x::InfiniteCardinal, ::Infinity) = x == ℵ₀
==(::Infinity, y::InfiniteCardinal) = ℵ₀ == y
==(::InfiniteCardinal{0}, y::RealInfinity) = ∞ == y
==(x::RealInfinity, ::InfiniteCardinal{0}) = x == ∞
==(::InfiniteCardinal, y::Real) = ∞ == y
==(x::Real, ::InfiniteCardinal) = x == ∞

@generated isless(::InfiniteCardinal{N}, ::InfiniteCardinal{M}) where {N,M} = :($(isless(N, M)))
isless(::InfiniteCardinal{0}, ::InfiniteCardinal{0}) = false
isless(x::Real, ::InfiniteCardinal{0}) = isfinite(x)
isless(x::Real, ::InfiniteCardinal) = true
isless(x::AbstractFloat, ::InfiniteCardinal{0}) = isfinite(x)
isless(x::AbstractFloat, ::InfiniteCardinal) = true
isless(::InfiniteCardinal, y::Real) = false
isless(x::InfiniteCardinal, y::AbstractFloat) = false

@generated <(::InfiniteCardinal{N}, ::InfiniteCardinal{M}) where {N,M} = :($(N < M))
@generated ≤(::InfiniteCardinal{N}, ::InfiniteCardinal{M}) where {N,M} = :($(N ≤ M))
@generated >(::InfiniteCardinal{N}, ::InfiniteCardinal{M}) where {N,M} = :($(N > M))
@generated ≥(::InfiniteCardinal{N}, ::InfiniteCardinal{M}) where {N,M} = :($(N ≥ M))

≤(::InfiniteCardinal{0}, ::InfiniteCardinal) = true
>(::InfiniteCardinal{0}, ::InfiniteCardinal) = false
≥(::InfiniteCardinal, ::InfiniteCardinal{0}) = true
<(::InfiniteCardinal, ::InfiniteCardinal{0}) = false

<(x::Real, ::InfiniteCardinal{0}) = x < ∞
<(x::Real, ::InfiniteCardinal) = true
≤(x::Real, ::InfiniteCardinal) = true
<(::InfiniteCardinal, x::Real) = false
≤(::InfiniteCardinal{0}, x::Real) = ∞ ≤ x
≤(::InfiniteCardinal, x::Real) = false
≤(::InfiniteCardinal{0}, x::RealInfinity) = ∞ ≤ x
≤(::InfiniteCardinal, x::RealInfinity) = false
>(::InfiniteCardinal{0}, y::Real) = ∞ > y
>(::InfiniteCardinal, ::Real) = true
≥(::InfiniteCardinal, ::Real) = true
>(::Real, ::InfiniteCardinal) = false
≥(x::Real, ::InfiniteCardinal{0}) = x ≥ ∞
≥(x::Real, ::InfiniteCardinal) = false

<(::Infinity, ::InfiniteCardinal{0}) = false
<(::Infinity, ::InfiniteCardinal) = true
≤(::Infinity, ::InfiniteCardinal) = true
<(::InfiniteCardinal, ::Infinity) = false
≤(::InfiniteCardinal{0}, ::Infinity) = true
≤(::InfiniteCardinal, ::Infinity) = false
>(::InfiniteCardinal{0}, ::Infinity) = false
>(::InfiniteCardinal, ::Infinity) = true
≥(::InfiniteCardinal, ::Infinity) = true
>(::Infinity, ::InfiniteCardinal) = false
≥(::Infinity, ::InfiniteCardinal{0}) = true
≥(::Infinity, ::InfiniteCardinal) = false

<(x::RealInfinity, ::InfiniteCardinal{0}) = x < ∞
<(x::RealInfinity, ::InfiniteCardinal) = true
≤(x::RealInfinity, ::InfiniteCardinal) = true
>(::InfiniteCardinal{0}, y::RealInfinity) = ∞ > y
>(::InfiniteCardinal, ::RealInfinity) = true
≥(::InfiniteCardinal, ::RealInfinity) = true

@generated min(::InfiniteCardinal{N}, ::InfiniteCardinal{M}) where {N,M} = :(InfiniteCardinal{$(min(N,M))}())
@generated max(::InfiniteCardinal{N}, ::InfiniteCardinal{M}) where {N,M} = :(InfiniteCardinal{$(max(N,M))}())
min(x::Real, ::InfiniteCardinal) = x
max(::Real, ℵ::InfiniteCardinal) = ℵ
min(::InfiniteCardinal, x::Real) = x
max(ℵ::InfiniteCardinal, ::Real) = ℵ

min(x::RealInfinity, ::InfiniteCardinal) = min(x, ∞)
max(::RealInfinity, ℵ::InfiniteCardinal) = ℵ
min(::InfiniteCardinal, y::RealInfinity) = min(∞, y)
max(ℵ::InfiniteCardinal, ::RealInfinity) = ℵ

min(::Infinity, ::InfiniteCardinal) = ∞
min(::InfiniteCardinal, ::Infinity) = ∞
max(::Infinity, ℵ::InfiniteCardinal) = ℵ
max(ℵ::InfiniteCardinal, ::Infinity) = ℵ

*(x::InfiniteCardinal) = x
*(::InfiniteCardinal{N}, ::InfiniteCardinal{N}) where N = InfiniteCardinal{N}()
*(::InfiniteCardinal{N}, ::Infinity) where N = InfiniteCardinal{N}()
*(::Infinity, ::InfiniteCardinal{N}) where N = InfiniteCardinal{N}()
function *(a::Integer, b::InfiniteCardinal)
    a > 0 || throw(ArgumentError("$a is non-positive"))
    b
end

*(a::Number, b::InfiniteCardinal) = a * ∞

*(a::InfiniteCardinal, b::Integer) = b*a
*(a::InfiniteCardinal, b::Number) = b*a

+(x::InfiniteCardinal) = x
+(x::InfiniteCardinal, y::InfiniteCardinal) = max(x,y)

+(::Integer, y::InfiniteCardinal) = y
+(x::InfiniteCardinal, ::Integer) = x
-(x::InfiniteCardinal, ::Integer) = x

-(::InfiniteCardinal{N}, ::InfiniteCardinal{N}) where N = NotANumber()

-(::InfiniteCardinal) = -∞
-(x::Integer, ::InfiniteCardinal) = x - ∞

for op in (:+, :-)
    @eval begin
        $op(x::Number, ::InfiniteCardinal) = $op(x, ∞)
        $op(::InfiniteCardinal, x::Number) = $op(∞, x)
        $op(x::RealInfinity, ::InfiniteCardinal) = $op(x, ∞)
        $op(::InfiniteCardinal, x::RealInfinity) = $op(∞, x)
        $op(::Infinity, ::InfiniteCardinal) = $op(∞, ∞)
        $op(::InfiniteCardinal, ::Infinity) = $op(∞, ∞)
    end
end

for OP in (:fld,:cld,:div)
    @eval begin
        $OP(x::InfiniteCardinal, ::Real) = x
        $OP(::InfiniteCardinal, ::InfiniteCardinal) = NotANumber()
    end
end

div(::T, ::InfiniteCardinal) where T<:Real = zero(T)
fld(x::T, ::InfiniteCardinal) where T<:Real = signbit(x) ? -one(T) : zero(T)
cld(x::T, ::InfiniteCardinal) where T<:Real = signbit(x) ? zero(T) : one(T)

mod(::InfiniteCardinal, ::InfiniteCardinal) = NotANumber()
mod(::InfiniteCardinal, ::Real) = NotANumber()
function mod(x::Real, ::InfiniteCardinal)
    x ≥ 0 || throw(ArgumentError("mod(x,∞) is unbounded for x < 0"))
    x
end

Base.to_index(::Union{Infinity,InfiniteCardinal{0}}) = ℵ₀
Base.to_shape(::Union{Infinity,InfiniteCardinal{0}}) = ℵ₀
Base.to_shape(dims::Tuple{Vararg{Union{Infinity, Integer, AbstractUnitRange}}}) = map(Base.to_shape, dims)

##
# Checked
##

Base.Checked.checked_sub(::Integer, ::InfiniteCardinal{0}) = -∞
Base.Checked.checked_sub(::InfiniteCardinal{0}, ::Integer) = ℵ₀
Base.Checked.checked_add(::Integer, ::InfiniteCardinal{0}) = ℵ₀
Base.Checked.checked_add(::InfiniteCardinal{0}, ::Integer) = ℵ₀

Base.Checked.checked_mul(x::Integer, ::InfiniteCardinal{0}) = sign(x)*∞
Base.Checked.checked_mul(::InfiniteCardinal{0}, x::Integer) = sign(x)*∞

##
# hash
##

Base.hash(::InfiniteCardinal{0}) = 0x775431eef01bca90 # made up
Base.hash(::InfiniteCardinal{1}) = 0x55437c69b794f8ce # made up

# avoid stack overflow
getindex(A::Array, i1::InfiniteCardinal{0}, I::Integer...) = throw(BoundsError(A, i1, I...))

Base._unsafe_getindex(::IndexStyle, A::AbstractArray, I::InfiniteCardinal{0}) = error("Overload getindex(::$(typeof(A)), ::InfiniteCardinal{0})")

# Avoid too-strict restrictions in SubArray
function getindex(V::SubArray{T,N}, i1::InfiniteCardinal{0}, I::Integer...) where {T,N}
    @boundscheck checkbounds(V, i1, I...)
    @inbounds r = V.parent[Base.reindex(V.indices, tuple(i1, I...))...]
    r
end


text:
module IntervalCensored

using Distributions
using Random
using SpecialFunctions
using HypergeometricFunctions
using ForwardDiff
using StaticArrays
using StatsFuns
import LinearAlgebra: dot
import Distributions: ccdf, cdf, logpdf, pdf, quantile,
       mean, rand, params, shape, scale

include("NonParametric.jl")
include("SimTools.jl")
include("calclp.jl")
include("./survdist/survdist.jl")
include("./survdist/GeneralizedGamma.jl")
include("./survdist/LogLogistic.jl")
include("MCEM.jl")

export calclp_ic, calclp_icrt, calclp_dic,
 make_ic, make_icrt, make_dic, make_dicrt,
 SurvIC, SurvICRT, SurvDIC, SurvDICRT,
 eqcdf, cdf, ccdf, pdf, logpdf, mean, shape, scale, params, rand,
 MCEMic, MCEMicrt, MCEMdic,
 GeneralizedGamma, LogLogistic

end

max_stars_repo_name: abikoushi/IntervalCensored

max_issues_repo_name: abikoushi/IntervalCensored

max_forks_repo_name: abikoushi/IntervalCensored


text:
using Nettle
using Test

include("hash_tests.jl")
include("hmac_tests.jl")
include("cipher_tests.jl")


text:
#==========================================================
This routine initializes an optimization problem to recover
a CNOT gate on a single qudit with 4 energy levels (and 2
guard states). The drift Hamiltonian in the rotating frame
is
              H_0 = - 0.5*ξ_a(a^†a^†aa),
where a is the annihilation operator for the qudit.
Here the control Hamiltonian includes the usual symmetric
and anti-symmetric terms
     H_{sym} = p(t)(a + a^†),    H_{asym} = q(t)(a - a^†)
which come from the rotating frame approximation and hence
we refer to these as "coupled" controls.
The problem parameters are:
                ω_a =  2π × 4.10595   Grad/s,
                ξ_a =  2π × 2(0.1099) Grad/s.
We use Bsplines with carrier waves with frequencies
0, ξ_a, 2ξ_a Grad/s.
==========================================================#
using LinearAlgebra
using Plots
using FFTW
using DelimitedFiles
using Printf
using Ipopt
using Random

Base.show(io::IO, f::Float64) = @printf(io, "%20.13e", f)

using Juqbox # quantum control module

Nosc = 1 # Number of coupled sub-systems = oscillators
N = 4 # Number of essential energy levels

Nguard = 2 # Number of extra guard levels
Ntot = N + Nguard # Total number of energy levels

T = 100.0 # Duration of gate.

# frequencies (in GHz, will be multiplied by 2*pi to get angular frequencies in the Hamiltonian matrix)
fa = 4.10336
xa = 0.2198
rot_freq = [fa] # Used to calculate the lab frame ctrl function

# setup drift Hamiltonian
number = Diagonal(collect(0:Ntot-1))

H0 = -0.5*(2*pi)*xa* (number*number - number) # xa is in GHz

# lowering matrix
amat = Bidiagonal(zeros(Ntot),sqrt.(collect(1:Ntot-1)),:U) # standard lowering matrix
adag = Array(transpose(amat));
Hsym_ops=[Array(amat + adag)]
Hanti_ops=[Array(amat - adag)]
H0 = Array(H0)

# Estimate time step
maxctrl = 0.001*2*pi * 8.5 #  9, 10.5, 12, 15 MHz

nsteps = calculate_timestep(T, H0, Hsym_ops, Hanti_ops, [maxctrl])
println("# time steps: ", nsteps)

Nfreq = 3 # number of carrier frequencies

Nctrl = length(Hsym_ops) # Here, Nctrl = 1
om = zeros(Nctrl,Nfreq)

om[1:Nctrl,2] .= -2.0*pi *xa       # Note negative sign
om[1:Nctrl,3] .= -2.0*pi* 2.0*xa

println("Carrier frequencies [GHz]: ", om[1,:]./(2*pi))

maxamp = zeros(Nfreq)

if Nfreq >= 3
    const_fact = 0.45
    maxamp[1] = maxctrl*const_fact
    maxamp[2:Nfreq] .= maxctrl*(1.0-const_fact)/(Nfreq-1) # split the remainder equally
else
    # same threshold for all frequencies
    maxamp .= maxctrl/Nfreq
end

maxpar = maximum(maxamp)

# Initial basis with guard levels
U0 = initial_cond([N], [Nguard])

# CNOT target
gate_cnot =  zeros(ComplexF64, N, N)
gate_cnot[1,1] = 1.0
gate_cnot[2,2] = 1.0
gate_cnot[3,4] = 1.0
gate_cnot[4,3] = 1.0

# Initial basis with guard levels
U0 = initial_cond([N], [Nguard])

utarget = U0 * gate_cnot

omega1 = Juqbox.setup_rotmatrices([N], [Nguard], [fa])

# Compute Ra*utarget
rot1 = Diagonal(exp.(im*omega1*T))

# target in the rotating frame
vtarget = rot1*utarget

params = Juqbox.objparams([N], [Nguard], T, nsteps, Uinit=U0, Utarget=vtarget, Cfreq=om, Rfreq=rot_freq,
                          Hconst=H0, Hsym_ops=Hsym_ops, Hanti_ops=Hanti_ops)

# initial parameter guess

# D1 smaller than 5 does not work
D1 = 10 # Number of B-spline coefficients per segment
nCoeff = 2*Nctrl*Nfreq*D1 # factor '2' is for sin/cos

# Random.seed!(12456)

startFromScratch = true # false
startFile="cnot-pcof-opt.dat"

# initial parameter guess
if startFromScratch
    pcof0 = maxpar*0.01 * rand(nCoeff)
    println("*** Starting from random pcof with amplitude ", maxpar*0.01)
else
    # use if you want to read the initial coefficients from file
    pcof0 = vec(readdlm(startFile))
    println("*** Starting from B-spline coefficients in file: ", startFile)
    nCoeff = length(pcof0)
    D1 = div(nCoeff, 2*Nctrl*Nfreq)  # number of B-spline coeff per control function
end

# min and max coefficient values
minCoeff, maxCoeff = Juqbox.assign_thresholds_freq(maxamp, Nctrl, Nfreq, D1)

samplerate = 32 # only used for plotting
casename = "cnot1" # base file name (used in optimize-once.jl)

maxIter = 75 # 0  # optional argument
lbfgsMax = 250 # optional argument
ipTol = 1e-5   # optional argument
acceptTol = ipTol # 1e-4 # acceptable tolerance
acceptIter = 15

println("*** Settings ***")
println("System Hamiltonian coefficients [GHz]: (fa, xa) =  ", fa, xa)
println("Total number of states, Ntot = ", Ntot, " Total number of guard states, Nguard = ", Nguard)
println("Using B-spline basis functions with carrier wave, # freq = ", Nfreq)
println("Number of coefficients per spline = ", D1, " Total number of parameters = ", nCoeff)
for q=1:Nfreq
    println("Carrier frequency: ", om[q]/(2*pi), " GHz, max parameter amplitude = ", 1000*maxamp[q]/(2*pi), " MHz")
end
println("Tikhonov coefficients: tik0 = ", params.tik0)

wa = Juqbox.Working_Arrays(params,nCoeff)
prob = Juqbox.setup_ipopt_problem(params, wa, nCoeff, minCoeff, maxCoeff, maxIter=maxIter, lbfgsMax=lbfgsMax, startFromScratch=startFromScratch)

println("Initial coefficient vector stored in 'pcof0'")


text:
#####
##### Interface for players
#####

"""
    AbstractPlayer{Game}

Abstract type for a player of `Game`.
"""
abstract type AbstractPlayer{Game} end

"""
    think(::AbstractPlayer, state, turn=nothing)

Return a probability distribution over actions as a `(actions, π)` pair.

The `turn` argument, if provided, indicates the number of actions that have
been played before by both players in the current game.
It is useful as during self-play, AlphaZero typically drops its temperature
parameter after a fixed number of turns.
"""
function think(::AbstractPlayer, state, turn=nothing)
  @unimplemented
end

"""
    reset_player!(::AbstractPlayer)

Reset the internal memory of a player (e.g. the MCTS tree).
The default implementation does nothing.
"""
function reset_player!(::AbstractPlayer)
  return
end

"""
    select_move(player::AbstractPlayer, state, turn=nothing)

Return a single action. A default implementation is provided that samples
an action according to the distribution computed by [`think`](@ref).
"""
function select_move(player::AbstractPlayer, state, turn=nothing)
  actions, π = think(player, state, turn)
  return actions[Util.rand_categorical(π)]
end

GameType(::AbstractPlayer{Game}) where Game = Game

#####
##### Random Player
#####

"""
    RandomPlayer{Game} <: AbstractPlayer{Game}

A player that picks actions uniformly at random.
"""
struct RandomPlayer{Game} <: AbstractPlayer{Game} end

function think(player::RandomPlayer, state, turn=nothing)
  actions = GI.available_actions(state)
  n = length(actions)
  π = ones(n) ./ length(actions)
  return actions, π
end

#####
##### Epsilon-greedy player
#####

"""
    EpsilonGreedyPlayer{Game, Player} <: AbstractPlayer{Game}

A wrapper on a player that makes it choose a random move
with a fixed ``ϵ`` probability.
"""
struct EpsilonGreedyPlayer{G, P} <: AbstractPlayer{G}
  player :: P
  ϵ :: Float64
  function EpsilonGreedyPlayer(p::AbstractPlayer{G}, ϵ) where G
    return new{G, typeof(p)}(p, ϵ)
  end
end

function think(p::EpsilonGreedyPlayer, state, turn=nothing)
  actions, π = think(p.player, state, turn)
  n = length(actions)
  η = ones(n) ./ n
  return actions, (1 - p.ϵ) * π + p.ϵ * η
end

function reset!(p::EpsilonGreedyPlayer)
  reset!(p.player)
end

#####
##### An MCTS-based player
#####

"""
    MctsPlayer{Game, MctsEnv} <: AbstractPlayer{Game}

A player that selects actions using MCTS.

# Constructors

    MctsPlayer(mcts::MCTS.Env; τ, niters, timeout=nothing)

Construct a player from an MCTS environment. When computing each move:

- if `timeout` is provided,
  MCTS simulations are executed for `timeout` seconds by groups of `niters`
- otherwise, `niters` MCTS simulations are run

The temperature parameter `τ` can be either a real number or a
[`StepSchedule`](@ref).

    MctsPlayer(oracle::MCTS.Oracle, params::MctsParams; timeout=nothing)

Construct an MCTS player from an oracle and an [`MctsParams`](@ref) structure.
If the oracle is a network, this constructor handles copying it, putting it
in test mode and copying it on the GPU (if necessary).
"""
struct MctsPlayer{G, M} <: AbstractPlayer{G}
  mcts :: M
  niters :: Int
  timeout :: Union{Float64, Nothing}
  τ :: StepSchedule{Float64} # Temperature
  function MctsPlayer(mcts::MCTS.Env{G}; τ, niters, timeout=nothing) where G
    @assert niters > 0
    @assert isnothing(timeout) || timeout > 0
    if isa(τ, Number)
      τ = StepSchedule(Float64(τ))
    else
      @assert isa(τ, StepSchedule)
    end
    new{G, typeof(mcts)}(mcts, niters, timeout, τ)
  end
end

# Alternative constructor
function MctsPlayer(
    oracle::MCTS.Oracle{G}, params::MctsParams; timeout=nothing) where G
  fill_batches = false
  if isa(oracle, AbstractNetwork)
    oracle = Network.copy(oracle, on_gpu=params.use_gpu, test_mode=true)
    params.use_gpu && (fill_batches = true)
  end
  mcts = MCTS.Env{G}(oracle,
    nworkers=params.num_workers,
    fill_batches=fill_batches,
    cpuct=params.cpuct,
    noise_ϵ=params.dirichlet_noise_ϵ,
    noise_α=params.dirichlet_noise_α)
  return MctsPlayer(mcts,
    niters=params.num_iters_per_turn,
    τ=params.temperature,
    timeout=timeout)
end

# MCTS with random oracle
function RandomMctsPlayer(::Type{G}, params::MctsParams) where G
  oracle = MCTS.RandomOracle{G}()
  mcts = MCTS.Env{G}(oracle,
    nworkers=1,
    cpuct=params.cpuct,
    noise_ϵ=params.dirichlet_noise_ϵ,
    noise_α=params.dirichlet_noise_α)
  return MctsPlayer(mcts,
    niters=params.num_iters_per_turn,
    τ=params.temperature)
end

function think(p::MctsPlayer, state, turn)
  if isnothing(p.timeout) # Fixed number of MCTS simulations
    MCTS.explore!(p.mcts, state, p.niters)
  else # Run simulations until timeout
    start = time()
    while time() - start < p.timeout
      MCTS.explore!(p.mcts, state, p.niters)
    end
  end
  return MCTS.policy(p.mcts, state, τ=p.τ[turn])
end

function reset_player!(player::MctsPlayer)
  MCTS.reset!(player.mcts)
end

#####
##### Network player
#####

"""
    NetworkPlayer{Game, Net} <: AbstractPlayer{Game}

A player that uses the policy output by a neural network directly,
instead of relying on MCTS.
"""
struct NetworkPlayer{G, N} <: AbstractPlayer{G}
  network :: N
  function NetworkPlayer(nn::AbstractNetwork{G}; use_gpu=true) where G
    nn = Network.copy(nn, on_gpu=use_gpu, test_mode=true)
    return new{G, typeof(nn)}(nn)
  end
end

function think(p::NetworkPlayer, state, turn)
  actions = GI.available_actions(state)
  board = GI.canonical_board(state)
  π, _ = MCTS.evaluate(p.network, board)
  return actions, π
end

#####
##### Players can play against each other
#####

"""
    play_game(white, black, memory=nothing)

Play a game between two [`AbstractPlayer`](@ref) and return the reward
obtained by `white`.

- If the `memory` argument is provided, samples are automatically collected
  from this game in the given [`MemoryBuffer`](@ref).
- If the `flip_probability` argument is set to ``p``, the board
  is _flipped_ randomly at every turn with probability ``p``,
  using [`GI.random_symmetric_state`](@ref).
"""
function play_game(
    white::AbstractPlayer{Game}, black::AbstractPlayer{Game},
    memory=nothing; flip_probability=0.) :: Float64 where Game
  state = Game()
  nturns = 0
  while true
    z = GI.white_reward(state)
    if !isnothing(z)
      isnothing(memory) || push_game!(memory, z, nturns)
      return z
    end
    if !iszero(flip_probability) && rand() < flip_probability
      state = GI.random_symmetric_state(state)
    end
    player = GI.white_playing(state) ? white : black
    actions, π = think(player, state, nturns)
    a = actions[Util.rand_categorical(π)]
    if !isnothing(memory)
      cboard = GI.canonical_board(state)
      push_sample!(memory, cboard, π, GI.white_playing(state), nturns)
    end
    GI.play!(state, a)
    nturns += 1
  end
end

#####
##### Evaluate two players against each other
#####

"""
    @enum ColorPolicy ALTERNATE_COLORS BASELINE_WHITE CONTENDER_WHITE

Policy for attributing colors in a duel between a baseline and a contender.
"""
@enum ColorPolicy ALTERNATE_COLORS BASELINE_WHITE CONTENDER_WHITE

"""
    pit(handler, contender, baseline, ngames)

Evaluate two `AbstractPlayer` against each other in a series of games.

# Arguments

  - `handler`: this function is called after each simulated
     game with two arguments: the game number `i` and the collected reward `z`
     for the contender player
  - `ngames`: number of games to play

# Optional keyword arguments
  - `reset_every`: if set, players are reset every `reset_every` games
  - `color_policy`: determines the [`ColorPolicy`](@ref),
     which is `ALTERNATE_COLORS` by default
  - `memory=nothing`: memory to use to record samples
  - `flip_probability=0.`: see [`play_game`](@ref)
"""
function pit(
    handler, contender::AbstractPlayer, baseline::AbstractPlayer, num_games;
    reset_every=nothing, color_policy=ALTERNATE_COLORS,
    memory=nothing, flip_probability=0.)
  baseline_white = (color_policy != CONTENDER_WHITE)
  zsum = 0.
  for i in 1:num_games
    white = baseline_white ? baseline : contender
    black = baseline_white ? contender : baseline
    z = play_game(white, black, memory, flip_probability=flip_probability)
    baseline_white && (z = -z)
    zsum += z
    handler(i, z)
    if !isnothing(reset_every) && (i % reset_every == 0 || i == num_games)
      reset_player!(baseline)
      reset_player!(contender)
    end
    if color_policy == ALTERNATE_COLORS
      baseline_white = !baseline_white
    end
  end
  return zsum / num_games
end

#####
##### Redudancy analysis
#####

# This type implements the interface of `MemoryBuffer` so that it can be
# passed as the `memory` argument of `play_game` to record games.
# Currently, this is only used to compute redundancy statistics
struct Recorder{Game, Board}
  boards :: Vector{Board}
  function Recorder{G}() where G
    B = GI.Board(G)
    return new{G, B}([])
  end
end

push_game!(r::Recorder, wr, gl) = nothing
push_sample!(r::Recorder, b, π, wp, t) = push!(r.boards, b)

function compute_redundancy(rec::Recorder{Game}) where Game
  initb = GI.board(Game())
  noninit = filter(!=(initb), rec.boards)
  unique = Set(noninit)
  return 1. - length(unique) / length(noninit)
end

#####
##### Simple utilities to play an interactive game
#####

"""
    Human{Game} <: AbstractPlayer{Game}

Human player that queries the standard input for actions.

Does not implement [`think`](@ref) but instead implements
[`select_move`](@ref) directly.
"""
struct Human{Game} <: AbstractPlayer{Game} end

struct Quit <: Exception end

function select_move(::Human, game, turn=nothing)
  a = nothing
  while isnothing(a) || a ∉ GI.available_actions(game)
    print("> ")
    str = readline()
    print("\n")
    isempty(str) && throw(Quit())
    a = GI.parse_action(game, str)
  end
  return a
end

"""
    interactive!(game, white, black)

Launch an interactive session for `game::AbstractGame` between players
`white` and `black`. Both players have type `AbstractPlayer` and one of them
is typically [`Human`](@ref).
"""
function interactive!(game, white, black)
  try
  GI.print_state(game)
  turn = 0
  while isnothing(GI.white_reward(game))
    player = GI.white_playing(game) ? white : black
    action = select_move(player, game, turn)
    GI.play!(game, action)
    GI.print_state(game)
    turn += 1
  end
  catch e
    isa(e, Quit) || rethrow(e)
    return
  end
end

interactive!(game::G) where G = interactive!(game, Human{G}(), Human{G}())


text:
# errors from LBL, remove post interpolation outliers

using Statistics

## Load necessary code

using Pkg
Pkg.activate(@__DIR__)
pkg"instantiate"
pkg"precompile"

using FileIO

include(joinpath(@__DIR__, "CommonUtils.jl"))
# include(joinpath(@__DIR__, "Plotting.jl"))

## load specific case

pargs = if !isdefined(Main, :parsed_args)
  parse_commandline()
else
  parsed_args
end

##

pargs["reportDir"] = if false && !haskey(pargs, "reportDir") && isdefined(Main, :fg)
  getLogPath(fg)
else
  # "/tmp/caesar/2020-02-23T01:43:32.222/fg_before_x291.tar.gz"
  # "/tmp/caesar/2020-02-24T04:22:23.282/fg_after_x1391.tar.gz"
  "/tmp/caesar/2020-02-24T04:19:39.943/fg_after_x1391.tar.gz"
end

# @show pargs["reportDir"]
# @show splitpath(pargs["reportDir"])

# load the factor graph needed
fg = if true
  println("going to load fg")
  fg = LightDFG{SolverParams}(params=SolverParams())
  @show pathElem = splitpath(pargs["reportDir"])
  @show getSolverParams(fg).logpath = joinpath(pathElem[1:end-1]...)
  loadDFG(pargs["reportDir"], Main, fg)
  fg
else
  fg
end

# random temp data
drt_data, TTm, XXm, YYm, XXf, YYf = loadResultsDRT(fg)

##

posesyms = ls(fg, r"x\d") |> sortDFG
filter!(x->isInitialized(fg, x), posesyms)
filter!(x->solverData(getVariable(fg, x), :lbl) != nothing, posesyms)
Xlbl = (x->(solverData(getVariable(fg, x), :lbl).val[1,1])).(posesyms)
Ylbl = (x->(solverData(getVariable(fg, x), :lbl).val[2,1])).(posesyms)

# pose PPE
Xppe = (x->(getVariablePPE(getVariable(fg, x)).suggested[1])).(posesyms)
Yppe = (x->(getVariablePPE(getVariable(fg, x)).suggested[2])).(posesyms)

# timestamp
ts = (x->getTimestamp(getVariable(fg, x))).(posesyms) .|> datetime2unix
T0 = ts[1]
ts .-= ts[1]
ts[end]

# filter on timestamps
# remove 420-450s
mask = 420 .< ts .< 450
imask = xor.(mask, true)

Xlbl[mask] .= Inf
Ylbl[mask] .= Inf

Gadfly.plot(x=ts, y=Xppe, Geom.line)
Gadfly.plot(x=ts, y=Yppe, Geom.line)

lblErr2 = (Xppe-Xlbl).^2 + (Yppe-Ylbl).^2
@show Statistics.mean(lblErr2[imask]) |> sqrt
@show cov(sqrt.(lblErr2[imask])) |> sqrt

##

Gadfly.plot(x=Xppe-Xlbl, Geom.histogram)
Gadfly.plot(x=ts, y=Xppe-Xlbl, Geom.line)
Gadfly.plot(x=ts, y=Yppe-Ylbl, Geom.line)

Gadfly.plot(x=ts, y=sqrt.(lblErr2), Geom.line)

findVariableNearTimestamp(fg,unix2datetime(T0+580))

## ranges nearby x580

ls(fg, :x581)

reportFactors(fg, Pose2Point2Range) #, [:x581l1f1; :x454l1f1])

ls(fg, :x454)

##  Plot LBL and belief only

drt_data, TTm, XXm, YYm, XXf, YYf = loadResultsDRT(fg)

mask = YYf .< -32
XXfm = XXf[mask]
YYfm = YYf[mask]
# only start does not match
TTmm = sum(mask)<length(TTm) ? TTm[end-sum(mask)+1:end] : TTm

# plot DRT
pl = Gadfly.plot(x=XXm[1:10:end], y=YYm[1:10:end], Geom.path) #, Theme(default_color=colorant"khaki"))

drtt = datetime2unix.(TTmm)
drtt .-= drtt[1]

pl = Gadfly.plot(
  Gadfly.layer(x=ts, y=Xppe, Geom.line, Theme(default_color=colorant"black")),
  Gadfly.layer(x=drtt, y=XXfm, Geom.line, Theme(default_color=colorant"red"))
)

#