text:
using Plots; gr()
using Reproduce
using FileIO
using Statistics
using ProgressMeter
using JLD2

# data_home = "../Experiment2_d_pi"
# data_home = "../data/OneDTMaze_RR_dpi"
# data_home = "../data/OneDTMaze_Control_dpi"
# data_home = "../OneDTMaze_RR_dpi"
# data_home = "../data/OneDTMaze_Control_dpi"
data_home = "../data/MC_Experiments"

include("./plot_utils.jl")
GPU = GeneralPlotUtils
LU = LabelUtils

experiment_folders = [data_home]
# folder_name = "oned_rr_dpi"
# folder_name = "oned_control_dpi"
folder_name = "MC_tmp"

# data_key = :oned_tmaze_dpi_error
# data_key = :ttmaze_direct_error
data_key = :mc_uniform_error
# data_key = :mc_starts_error

ic = ItemCollection(joinpath(experiment_folders[1], "data"));
ic = search(ic, Dict("demon_learner" => "SR"))

algo_divisor_keys = ["behaviour_learner", "demon_learner"]
sweep_params = ["demon_eta","behaviour_eta","exploration_param"]
# algo_divisor_keys = ["behaviour_learner", "demon_learner", "demon_opt", "demon_update", "intrinsic_reward"]
# sweep_params = ["demon_alpha_init", "demon_eta"]
#
algo_specs_full = GPU.split_algo(ic, algo_divisor_keys)

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

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

algo_specs = algo_specs_full[valid_algos_ind]
algo_ics = all_algos_ics[valid_algos_ind]
println()
@show algo_specs

best_per_algo_ics = []
for (i,algo_ic) in enumerate(algo_ics)
#     push!(best_per_algo_ics, GPU.get_best(algo_ic,sweep_params, data_key))
    push!(best_per_algo_ics, GPU.get_best_final_perf(algo_ic,sweep_params, data_key, 0.1))
    # push!(best_per_algo_ics, GPU.get_most_episodes(algo_ic,sweep_params, :episode_length))
end
@show length.(best_per_algo_ics), "better be equal to num seeds"

# color = LU.get_colour(algo_ics[1])
# style = LU.get_linestyle(algo_ics[1])

plot_params = [LU.get_params(best_per_algo_ics[i]) for i in 1:length(best_per_algo_ics)]
# plot(rand(10); plot_params...)
#
num_runs = length.(best_per_algo_ics)[1]
@show num_runs

# data = [GPU.smooth(GPU.get_stats(GPU.load_results(ic,data_key))[1],10) for ic in best_per_algo_ics]
# std = [GPU.smooth(GPU.get_stats(GPU.load_results(ic,data_key))[2],10) for ic in best_per_algo_ics]

data = [GPU.smooth(GPU.get_stats(GPU.load_results(ic,data_key))[1],10) for ic in best_per_algo_ics]
std = [GPU.smooth(GPU.get_stats(GPU.load_results(ic,data_key))[2],10) for ic in best_per_algo_ics]
#
# [GPU.print_params(algo, algo_divisor_keys, sweep_params) for algo in best_per_algo_ics]
# label_keys = cat(algo_divisor_keys,sweep_params, dims = 1)
# labels = [GPU.get_label(algo, label_keys) for algo in best_per_algo_ics]
# labels = cat(labels..., dims=2)
# string(data_key)
#
# ####
# Generate Average demon RMSE
ylabel = "RMSE"
# title = "SR Demons & Step Size Adaptation vs More Naive Approaches"
title = "RMSE"
xlabel = "Steps"
step_increment=1000
num_samples = length(data[1])
xticks=collect(1:step_increment:num_samples*step_increment)
p = plot(ylabel=ylabel, grid=true, title=title, xlabel=xlabel)
for i in 1:length(best_per_algo_ics)
    l = GPU.get_label(best_per_algo_ics[i],cat(algo_divisor_keys, sweep_params, dims = 1))
    plot!(p, xticks, data[i], ribbon = std[i]/sqrt(num_runs), legend=:topright, label = l)
end
savefig("./plots/$(folder_name)/RMSE_$(string(data_key)).png")

# ####
# # Generate graph for RMSE per demon
#
data_per_gvf = [GPU.get_stats(GPU.load_results(ic,data_key), per_gvf = true)[1] for ic in best_per_algo_ics]
std_per_gvf = [GPU.get_stats(GPU.load_results(ic,data_key), per_gvf = true)[2] for ic in best_per_algo_ics]
# num_gvfs = 4
num_gvfs = 2

ps = []
# gvf_labels = ["Distractor" "Constant" "Drifter" "Constant"]
gvf_labels = ["Wall", "Goal"]

for gvf_ind in 1:num_gvfs
    p = plot()
    for algo_ind in 1:length(data_per_gvf)
        smooth_gvf = GPU.smooth(data_per_gvf[algo_ind][gvf_ind,:],5)
        label = GPU.get_label(best_per_algo_ics[algo_ind],algo_divisor_keys)

        plot!(p,smooth_gvf, palette=:tab10, ribbon = std_per_gvf[algo_ind][gvf_ind,:] / sqrt(num_runs), size = (500,500),label = label)
        plot!(xlabel="Steps", ylabel = "RMSE",  title = string(gvf_labels[gvf_ind], ""))
    end
    push!(ps,p)
end
plot(ps..., size = (1200,1200))
savefig("./plots/$(folder_name)/RMSE_per_Demon_$(string(data_key)).png")
#
# # Generate Goal visitation
episode_lengths = [GPU.load_results(ic,:episode_length, return_type = "array") for ic in best_per_algo_ics]
# visit_orders = [GPU.load_results(ic,:visit_order, return_type = "array") for ic in best_per_algo_ics]
# ps = []
# for ind in 1:length(episode_lengths)
#     num_episodes = [length(e) for e in episode_lengths[ind]]
#     p = histogram(num_episodes; nbins=15, xlim=(0,3000))
#
#     title!(labels[ind])
#     push!(ps,p)
#     xlabel!("Number of Episodes")
#     ylabel!("Number of runs")
# end
#
# plot(ps..., size = (1000,1000))
# savefig("./plots/$(folder_name)/RunSuccessRate_ALL_$(string(data_key)).png")
#
# #### Generate
max_lengths = [GPU.get_min_length(arrs)-1 for arrs in episode_lengths]
@show max_lengths
@show episode_lengths
p = plot(xlabel = "Episode Count", ylabel = "Step Length")
for (ind,episode_lengths) in enumerate(episode_lengths)
    l = hcat([epi[1:max_lengths[ind]] for epi in episode_lengths]...)
    label = GPU.get_label(best_per_algo_ics[ind],algo_divisor_keys)
    plot!(p,mean(l,dims=2), label = label)
end
title!("Step Length for Best Performing Algos Last 10 %")
# p = plot(mean(episode_lengths,dims=2), xlabel="Episode Count", ylabel="Step Length", hline=5.5)
display(p)
savefig("./plots/$(folder_name)/step_length.png")
#
# ####
# ps = []
# for (algo_ind,algo) in enumerate(episode_lengths)
#     p = plot()
#     num_runs = length(algo)
#     num_steps = cumsum(algo[1])[end]
#     d = zeros(num_runs, num_steps)
#     for (i,run) in enumerate(algo)
#         cur_index = 1
#         for (ep_num,ep) in enumerate(run)
#             d[i,cur_index:cur_index + ep-1] .= ep
#             cur_index += ep
#         end
#             plot!(d[i,:],color=:black, linealpha=0.4, label="")
#     end
#     plot!(mean(d,dims=1)', label = labels[algo_ind], yaxis=:log, size=(1000,1000), legend=:bottomleft)
#     xlabel!("Step")
#     ylabel!("Ave Ep Length at Step")
#     push!(ps, p)
#
# end
# plot(ps...)
# savefig("./plots/$(folder_name)/ep_length_per_step_$(string(data_key)).png")
#
# #####
# # Goal Visitation
# ps = []
# for i in 1:length(visit_orders)
#     @show max_lengths[i]
#     @show length(visit_orders[i])
#     visit_perc = GPU.goal_visits_per_episode(visit_orders[i], max_lengths[i])
#     @show size(visit_perc)
#     tmp = [GPU.smooth(visit_perc[i,:],20) for i in 1:4]
#     p = plot(tmp, labels = gvf_labels, xlabel="Episode Count", ylabel="Fraction of Goal Visits", ylim=(0.0,1.0), title = string(labels[i], " Last 10%"), xaxis=:log, size = (800,1000), legend=:topleft)
#     push!(ps,p)
# end
# println()
# @show length(visit_orders)
# plot(ps..., layout=(Int(length(best_per_algo_ics)/2),2))
# savefig("./plots/$(folder_name)/goal_visits_$(string(data_key)).png")
#
# # Median Goal Visitation
# tmp_lengths = [length.(visits) for visits in visit_orders]
# med_visit_order_lengths = Int.(floor.([median(l) for l in tmp_lengths]))
# med_visit_mask = [tmp_lengths[i] .> med_visit_order_lengths[i] for i in 1:length(visit_orders)]
# ps = []
#
# for i in 1:length(visit_orders)
#     visit_perc = GPU.goal_visits_per_episode(visit_orders[i][med_visit_mask[i]], med_visit_order_lengths[i])
#     tmp = [GPU.smooth(visit_perc[i,:],100) for i in 1:4]
#     p = plot(tmp, labels = gvf_labels, xlabel="Episode Count", ylabel="Fraction of Goal Visits", ylim=(0.0,1.0), title = string("MEDIAN",labels[i], " Last 10%"), size = (800,1000), legend=:topleft)
#     push!(ps,p)
# end
# println()
# @show length(visit_orders)
# display(plot(ps..., layout=(Int(length(best_per_algo_ics)/2),2)))
# savefig("./plots/$(folder_name)/goal_visits_MEDIAN_$(string(data_key)).png")

max_stars_repo_name: MatthewMcLeod/curiosity

max_issues_repo_name: MatthewMcLeod/curiosity

max_forks_repo_name: MatthewMcLeod/curiosity


text:
module InterfacesTests

include("mpiexec.jl")
run_mpi_driver(procs=4,file="driver_interfaces.jl")

end # module

ext: jl


text:
const OpenClosed{T} = IntervalSets.TypedEndpointsInterval{:open,:closed, T}
const ClosedOpen{T} = IntervalSets.TypedEndpointsInterval{:closed,:open, T}
const OpenOpen{T} = IntervalSets.TypedEndpointsInterval{:open,:open, T}
const ClosedClosed{T} = IntervalSets.TypedEndpointsInterval{:closed,:closed, T}

_isless(a,b) = Base.isless(a,b)
_isless(a::AbstractInterval, b) = _isless(leftendpoint(a), b)

function Base.isless(a::Interval{:closed,:open}, b::IntervalSets.Interval{:closed,:open})
    Base.isless(leftendpoint(a), leftendpoint(b))
end


text:
# MODULARPLOT.jl
# Plots modular bipartite network configuration using row and column module labels. Row and column numbers indicate the ordering of rows and columns within the original network.
# Plotting performed using the PyPlot library.
# Author :  Stephen Beckett ( https://github.com/sjbeckett/weighted-modularity-LPAwbPLUS )
# MIT License

function MODULARPLOT(MATRIX,OUT)
    #MATRIX:  original bipartite network
    #OUT: output from modularity algorithm containing row and column module labels
    print("Attempting to import PyPlot library...")
    require("PyPlot")
    print("Plotting...")
    #1. Sort original network, so that rows are positioned according to modules
    ROW_IX = sortperm(vec(OUT[1]));
    COL_IX = sortperm(vec(OUT[2]));
    ROWS = OUT[1][ROW_IX];
    COLS = OUT[2][COL_IX];
    MODS = unique(ROWS);

    #2. Plot this network
    MATRIXPLOT(MATRIX[ROW_IX,COL_IX], COL_IX, ROW_IX)

    #Add the modular boundaries using output information - drawing from top left to bottom right
    xcurr=0.5
    ycurr=length(ROWS)+0.5

    for aa = 1:length(MODS)
        Rsize = sum(ROWS.==MODS[aa]) #Number of rows in this module
        Csize = sum(COLS.==MODS[aa]) #Number of columns in this module
        PyPlot.fill_betweenx([ycurr-0.5,ycurr-Rsize-0.5],xcurr-0.5,xcurr+Csize-0.5,color="none",edgecolor="red")
        xcurr = xcurr+Csize
        ycurr = ycurr-Rsize
    end

end

###############################################

function MATRIXPLOT(INMATRIX,col_ix,row_ix)
	BACKGROUNDCOLOR = (173/255,235/255,255/255);
	offset = 0.3;

	xs = length(col_ix)
	ys = length(row_ix)

	#initial background
	PyPlot.pcolor(INMATRIX)
	PyPlot.xticks((1:xs)-0.5,col_ix)
	PyPlot.yticks((1:ys)-0.5,reverse(row_ix))
	PyPlot.fill_betweenx([0:ys],0,xs,color=BACKGROUNDCOLOR)

	#fill in
	minR = minimum(INMATRIX)
	RANGE = maximum(INMATRIX)-minR
	scalefactor = 1/RANGE

	for b = 1:xs
		for a = 1:ys
			if INMATRIX[1+ys-a,b]!=0
				COL = 1-(INMATRIX[1+ys-a,b]-minR)*scalefactor
				PyPlot.fill_betweenx([a-offset-0.5,a+offset-0.5],b-offset-0.5,b+offset-0.5,color=(COL,COL,COL))
			end
		end
	end
end

max_stars_repo_name: sjbeckett/weighted-modularity-LPAwbPLUS

max_issues_repo_name: sjbeckett/weighted-modularity-LPAwbPLUS

max_forks_repo_name: sjbeckett/weighted-modularity-LPAwbPLUS


text:
"""
    zmfillw()

Fill weight tables for zero-mass structure function
calculations.
"""
function zmfillw()

    nwords = Ref{Int32}()

    @ccall zmfillw_(nwords::Ref{Int32})::Nothing

    nwords[]
end

"""
    zmwords()

Check the number of words available in the ZMSTF
workspace and the number of words used.
"""
function zmwords()

    ntotal = Ref{Int32}()
    nused = Ref{Int32}()

    @ccall zmwords_(ntotal::Ref{Int32}, nused::Ref{Int32})::Nothing

    ntotal[], nused[]
end

"""
    zmstfun()

Calculate a structure function from a linear combination
of parton densities.

# Arguments
- `istf::Integer`: structure function index
where (1,2,3,4) = (F_L, F_2, xF_3, f_L^').
- `def::Array{Float64}`: coeffs of the quark linear combination
for which the structure function is to be calculated.
- `x::Array{Float64}`: list of `x` values.
- `Q2::Array{Float64}`: list of `Q2` values.
- `n::Integer`: number of items in `x`, `Q2` and `f`.
- `ichk::Integer`: flag for grid boundary checks. See QCDNUM docs.

# Returns
- `f::Array{Float64}`: list of structure functions.
"""
function zmstfun(istf::Integer, def::Array{Float64}, x::Array{Float64},
                 Q2::Array{Float64}, n::Integer, ichk::Integer)

    istf = Ref{Int32}(istf)
    n = Ref{Int32}(n)
    ichk = Ref{Int32}(ichk)

    f = Array{Float64}(undef, n[])

    @ccall zmstfun_(istf::Ref{Int32}, def::Ref{Float64}, x::Ref{Float64},
                    Q2::Ref{Float64}, f::Ref{Float64}, n::Ref{Int32},
                    ichk::Ref{Int32})::Nothing

    f
end


text:
## -------------------------------------------------------------------
# flx correlations
let

    MEname(c) = "\$\\textbf{ME}^{\\langle \\textbf{$c} \\rangle} \\textbf{ (abs. flux)} \$"

    me_method = :ME_MAX_POL
    fba_method = :FBA_Z_FIX_MIN_COST
    me_label = MEname(1)
    fba_label = "\$ \\textbf{FBA} \\textbf{ (abs. flux)} \$"

    for (get_flxs, flxs) in [
        ((DAT) -> DAT[:EXCH_FLX_IDERS], "all"),
        ((DAT) -> ["GLC", "AC"], "common"),
    ]

        fontsize = 13
        p = plot(;
            title = "exchanges",
            xlabel = fba_label, ylabel = me_label,
            thickness_scaling = 1.3, dpi = 1000,
            xguidefontsize = fontsize, yguidefontsize = fontsize
        )

        all_vals = []
        for (iJR, Data) in [
            (ChK.iJR904, ChK.KayserData),
            (ChN.iJR904, ChN.NanchenData),
            (ChF.iJR904, ChF.FolsomData),
            # (ChH.iJR904, ChH.HeerdenData),
        ]
            src = nameof(Data)
            marker = (6, source_markers[Data])

            datfile = iJR.procdir("dat.bson")
            DAT = UJL.load_data(datfile; verbose = false)

            # EXCH_FLX_IDERS = DAT[:EXCH_FLX_IDERS]
            # EXCH_FLX_IDERS = ["GLC", "AC"]
            EXCH_FLX_IDERS = get_flxs(DAT)
            EXPS = DAT[:EXPS]

            tot_ps = Plots.Plot[]

            # total corr
            let
                ep_vals = DAT[me_method, :flx, EXCH_FLX_IDERS, EXPS] .|> abs
                ep_errs = DAT[me_method, :err, EXCH_FLX_IDERS, EXPS] .|> abs

                fba_vals = DAT[fba_method, :flx, EXCH_FLX_IDERS, EXPS] .|> abs

                color = [ider_colors[ider] for ider in EXCH_FLX_IDERS, exp in EXPS]

                # ep corr
                scatter!(p, fba_vals, ep_vals; yerr = ep_errs,
                    marker, label = "", color, alpha = 0.7
                )
                push!(all_vals, [ep_vals; fba_vals])
            end
        end

        sort!(all_vals)
        plot!(p, all_vals, all_vals; alpha = 0.7,
            ls = :solid, color = :black, label = ""
        )
        mysavefig(p, "fba_vs_me_exch_corr"; flxs)
    end
end


text:
#!/usr/bin/env julia
using CSV
using StatsPlots

df1 = CSV.read("julia_times.csv")
df2 = CSV.read("python_times.csv")

# rename systems
df1[:system] = ["jl_" * s for s in df1[:system]]
df2[:system] = ["py_" * s for s in df2[:system]]
df = vcat(df1, df2)

# Plot
@df df plot(
    :n,
    :time,
    group = :system,
    markershape = :circle,
    markersize = 4,
    markerstrokealpha = 0,
    linewidth = 1,
    alpha = 0.8,
    xlabel = "Number of coordinates",
    ylabel = "Time (s)",
    xscale = :log10,
    yscale = :log10,
    legend = :bottomright,
)
savefig("bench.png")


text:
module ClusteredLowRankSolver

using AbstractAlgebra, IterTools, LinearAlgebra
# using ApproximateFekete
using Printf, Arblib, BlockDiagonals
using KrylovKit

export LowRankMat, LowRankMatPol, Block, Constraint, Objective,	 LowRankPolProblem, ClusteredLowRankSDP, solvesdp, approximatefekete,SampledMPolyElem,optimal
export solvesdp, convert_to_prec, SolverFailure

import LinearAlgebra: dot, transpose
import Base: ==

import AbstractAlgebra: evaluate

include("interface.jl")
include("tools.jl")
include("threadinginfo.jl")
include("solver.jl")
include("approximate_fekete.jl") #I think this is better than having another auxiliary package

end


text:
function _precompile_()
    ccall(:jl_generating_output, Cint, ()) == 1 || return nothing
    Base.precompile(Tuple{typeof(servedocs_callback!),SimpleWatcher,String,String})   # time: 0.31880945
    isdefined(LiveServer, Symbol("#12#17")) && Base.precompile(Tuple{getfield(LiveServer, Symbol("#12#17")),HTTP.Streams.Stream{HTTP.Messages.Request, HTTP.ConnectionPool.Transaction{TCPSocket}}})   # time: 0.098543696
    Base.precompile(Tuple{typeof(get_fs_path),String})   # time: 0.087591365
    Base.precompile(Tuple{Core.kwftype(typeof(serve)),NamedTuple{(:port,), Tuple{Int64}},typeof(serve)})   # time: 0.07307711
    isdefined(LiveServer, Symbol("#11#16")) && Base.precompile(Tuple{getfield(LiveServer, Symbol("#11#16"))})   # time: 0.06859666
    isdefined(LiveServer, Symbol("#2#3")) && Base.precompile(Tuple{getfield(LiveServer, Symbol("#2#3"))})   # time: 0.030429186
    Base.precompile(Tuple{typeof(ws_tracker),HTTP.WebSockets.WebSocket,AbstractString})   # time: 0.03029817
    isdefined(LiveServer, Symbol("#6#7")) && Base.precompile(Tuple{getfield(LiveServer, Symbol("#6#7")),HTTP.WebSockets.WebSocket{IOBuffer}})   # time: 0.018438155
    Base.precompile(Tuple{typeof(watch_file!),SimpleWatcher,String})   # time: 0.0108271
    Base.precompile(Tuple{typeof(ws_upgrade),HTTP.Streams.Stream{HTTP.Messages.Request, IOBuffer}})   # time: 0.010094539
    Base.precompile(Tuple{typeof(set_callback!),SimpleWatcher,Function})   # time: 0.003649323
    Base.precompile(Tuple{typeof(append_slash),String})   # time: 0.003343235
    Base.precompile(Tuple{typeof(stop),SimpleWatcher})   # time: 0.002390376
    Base.precompile(Tuple{Type{SimpleWatcher}})   # time: 0.002108754
    Base.precompile(Tuple{typeof(start),SimpleWatcher})   # time: 0.001990792
    Base.precompile(Tuple{Core.kwftype(typeof(Type)),NamedTuple{(:sleeptime,), Tuple{Float64}},Type{SimpleWatcher},typeof(identity)})   # time: 0.00174082
    Base.precompile(Tuple{typeof(example)})   # time: 0.001634037
    Base.precompile(Tuple{typeof(has_changed),WatchedFile{String}})   # time: 0.001403002
    Base.precompile(Tuple{typeof(ws_tracker),HTTP.WebSockets.WebSocket{IOBuffer},String})   # time: 0.001390946
end


text:
v = [1, 2, 3, 4]
dv = DataArray(v, falses(size(v)))

m = [1 2; 3 4]
dm = DataArray(m, falses(size(m)))

t = Array(Int, 2, 2, 2)
t[1:2, 1:2, 1:2] = 1
dt = DataArray(t, falses(size(t)))

dv = DataArray(v)
dv = DataArray(v, [false, false, false, false])

dv = DataArray(dv)
dm = DataArray(dm)
dt = DataArray(dt)

dv = DataArray(Int, 2)
dm = DataArray(Int, 2, 2)
dt = DataArray(Int, 2, 2, 2)

similar(dv)
similar(dm)
similar(dt)

similar(dv, 2)
similar(dm, 2, 2)
similar(dt, 2, 2, 2)


text:
# ---
# jupyter:
#   jupytext:
#     formats: ipynb,jl:hydrogen
#     text_representation:
#       extension: .jl
#       format_name: hydrogen
#       format_version: '1.3'
#       jupytext_version: 1.11.2
#   kernelspec:
#     display_name: Julia 1.8.0-DEV
#     language: julia
#     name: julia-1.8
# ---

# %%
using StatsPlots
using Distributions

# %%
n = 10^3
bin = Binomial(n, 5/n)
k = 0:15
plot(k, pdf.(bin, k))
plot!(k, pdf.(Normal(mean(bin), std(bin)), k))
plot!(k, pdf.(Poisson(5), k))

# %%

max_stars_repo_path: 0013/Poisson approx of Binomial.jl

max_stars_repo_name: genkuroki/public

max_issues_repo_path: 0013/Poisson approx of Binomial.jl

max_issues_repo_name: genkuroki/public

max_forks_repo_path: 0013/Poisson approx of Binomial.jl

max_forks_repo_name: genkuroki/public


text:
##
using Random, LinearAlgebra
Random.seed!(0)
X1 = randn(4,10)
DX = zeros(size(X1,2),size(X1,2))
for i=1:size(X1,2)
  for j=1:size(X1,2)
    DX[i,j] = norm(X1[:,i]-X1[:,j])
  end
end
##
function nnsimple(X,k)
  m,n = size(X)
  idxs = zeros(Int,k,n)
  dsts = zeros(Float64,k,n)
  d = zeros(n)
  ix = ones(Int,n)
  subset = 1:k
  @inbounds for i=1:n
    fill!(d,0)
    for j=1:n
      @simd for k=1:m
        d[j] += (X[k,i] - X[k,j])^2
      end
      d[j] = sqrt(d[j])
    end
    partialsortperm!(ix,d,subset);
    for j=1:k
      idxs[j,i] = ix[j]
      dsts[j,i] = d[ix[j]]
    end
  end
  return idxs,dsts
end
begin
  idxs,dsts = nnsimple(X1,10)
  for i=1:size(X1,2)
    @assert dsts[:,i] ≈ sort(DX[:,i])
  end
end
using SparseArrays
function nngraph(idxs::Matrix{T},dsts::Matrix{F}) where {T <: Integer, F <: Real}
  n = size(idxs, 2)
  ei = zeros(Int,0)
  ej = copy(ei)
  ev = zeros(Float64,0)
  for i=1:n
    for k=1:size(idxs,1)
      if i != idxs[k,i]
        push!(ei, i)
        push!(ej, idxs[k,i])
        push!(ev, dsts[k,i])

        push!(ei, idxs[k,i])
        push!(ej, i)
        push!(ev, dsts[k,i])
      end
    end
  end
  A = sparse(ei,ej,ev,n,n,max)
end

max_stars_repo_name: empet/more-graphs

max_issues_repo_name: empet/more-graphs

max_forks_repo_name: empet/more-graphs


text:
@doc raw"""
    DAR(y::Vector{Fl}, lags::Int) where Fl

A Dynamic Autorregressive model is defined by:
```math
\begin{gather*}
    \begin{aligned}
        y_{t} &= \phi_0 + \sum_{i=1}^{lags} \phi_{i, t} y_{t - i} + \varepsilon_{t}  \quad \varepsilon_{t} \sim \mathcal{N}(0, \sigma^2_{\varepsilon})\\
        \phi_{i, t+1} &= \phi_{i, t} + \eta{i, t}  \quad \eta{i, t} \sim \mathcal{N}(0, \sigma^2_{i, \eta})\\
    \end{aligned}
\end{gather*}
```

!!!! note System matrices
    When building the system matrices we get rid of the first `lags` observations
    in order to build all system matrices respecting the correspondent timestamps

!!! warning Forecasting
    The dynamic autorregressive model does not have the [`forecast`](@ref) method implemented yet.
    If you wish to perform forecasts with this model please open an issue.

# Example
```jldoctest
julia> model = DAR(randn(100), 2)
DAR
```
"""
mutable struct DAR <: StateSpaceModel
    lags::Int
    first_observations::Vector{<:Real}
    hyperparameters::HyperParameters
    system::LinearUnivariateTimeVariant
    results::Results

    function DAR(y::Vector{Fl}, lags::Int) where Fl

        X = lagmat(y, lags)
        num_observations = size(X, 1)
        first_observations = y[1:lags]

        # Define system matrices
        Z = [X[t, :] for t in 1:num_observations]
        T = [Matrix{Fl}(I, lags, lags) for _ in 1:num_observations]
        R = [Matrix{Fl}(I, lags, lags) for _ in 1:num_observations]
        d = [zero(Fl) for _ in 1:num_observations]
        c = [zeros(lags) for _ in 1:num_observations]
        H = [one(Fl) for _ in 1:num_observations]
        Q = [Matrix{Fl}(I, lags, lags) for _ in 1:num_observations]

        system = LinearUnivariateTimeVariant{Fl}(y[lags+1:end], Z, T, R, d, c, H, Q)

        # Define hyperparameters names
        names = [["sigma2_ε", "intercept"];["sigma2_ar_$i" for i in 1:lags]]
        hyperparameters = HyperParameters{Fl}(names)

        return new(lags, first_observations, hyperparameters, system, Results{Fl}())
    end
end

# Obligatory methods
function default_filter(model::DAR)
    Fl = typeof_model_elements(model)
    a1 = zeros(Fl, num_states(model))
    # P1 is identity on sigma \eta and 0 otherwise.
    P1 = Matrix{Fl}(I, num_states(model), num_states(model))
    steadystate_tol = Fl(1e-5)
    return UnivariateKalmanFilter(a1, P1, 0, steadystate_tol)
end

function get_sigma_ar_name(model::DAR, i::Int)
    return model.hyperparameters.names[i + 2]
end

function initial_hyperparameters!(model::DAR)
    Fl = typeof_model_elements(model)
    observed_variance = var(model.system.y[findall(!isnan, model.system.y)])
    observed_mean = mean(model.system.y[findall(!isnan, model.system.y)])
    initial_hyperparameters = Dict{String,Fl}(
        "sigma2_ε" => observed_variance, "intercept" => observed_mean
    )
    for i in axes(model.system.T[1], 2)
        initial_hyperparameters[get_sigma_ar_name(model, i)] = one(Fl)
    end
    set_initial_hyperparameters!(model, initial_hyperparameters)
    return model
end

function constrain_hyperparameters!(model::DAR)
    for i in axes(model.system.T[1], 2)
        constrain_variance!(model, get_sigma_ar_name(model, i))
    end
    constrain_variance!(model, "sigma2_ε")
    constrain_identity!(model, "intercept")
    return model
end

function unconstrain_hyperparameters!(model::DAR)
    for i in axes(model.system.T[1], 2)
        unconstrain_variance!(model, get_sigma_ar_name(model, i))
    end
    unconstrain_variance!(model, "sigma2_ε")
    unconstrain_identity!(model, "intercept")
    return model
end

function fill_model_system!(model::DAR)
    H = get_constrained_value(model, "sigma2_ε")
    fill_H_in_time(model, H)
    intercept = get_constrained_value(model, "intercept")
    for t in 1:length(model.system.Q)
        fill!(model.system.d, intercept)
        for i in axes(model.system.Q[1], 1)
            model.system.Q[t][i, i] = get_constrained_value(model, get_sigma_ar_name(model, i))
        end
    end
    return model
end

function fill_H_in_time(model::DAR, H::Fl) where Fl
    return fill_system_matrice_with_value_in_time(model.system.H, H)
end

function reinstantiate(model::DAR, y::Vector{Fl}) where Fl
    return DAR(y, model.lags)
end

# This model needs a custom made forecasting routine
function forecast(
    model::DAR, steps_ahead::Int; filter::KalmanFilter=default_filter(model)
)
    return error("Forecasting is currently nor implemented for the DAR model.")
end


text:
function lyapdiag!(A::Vector,B::Vector,F)
    for i=1:size(F,1),j=1:size(F,2)
        F[i,j] *= 1/(A[i]+B[j])
    end

    F
end

# Solve P*Y*R' + S*Y*T' = F
# where P,R,S and T are upper triangular
function lyapuptriang(P,R,S,T,F::AbstractArray{N}) where N
    m=size(P,2); n = size(T,2)
    Y=Matrix{N}(m,n)
    PY = Matrix{N}(m,n); SY = Matrix{N}(m,n)

    k=n
    while k > 1
        if T[k,k-1] == 0 && R[k,k-1] == 0
            rhs = F[:,k]
            if k < n
                for j = k+1:n, l=1:m
                    rhs[l] -= R[k,j]*PY[l,j] + T[k,j]*SY[l,j]
                end
            end

            Y[:,k] = (R[k,k]*P + T[k,k]*S)\rhs
            PY[:,k] = P*Y[:,k];SY[:,k] = S*Y[:,k]

            k -= 1
        else
            rhs1 = F[:,k-1];rhs2=F[:,k]
            for j = k+1:n, l=1:m
                rhs1[l] -= R[k-1,j]*PY[l,j] + T[k-1,j]*SY[l,j]
                rhs2[l] -= R[k,j]*PY[l,j] + T[k,j]*SY[l,j]
            end

            SM = Matrix{N}(2m,2m)

            for i=1:m,j=1:m
                SM[i,j] = R[k-1,k-1]*P[i,j] + T[k-1,k-1]*S[i,j]
                SM[i,j+m] = R[k-1,k]*P[i,j] + T[k-1,k]*S[i,j]
                SM[i+m,j] = R[k,k-1]*P[i,j] + T[k,k-1]*S[i,j]
                SM[i+m,j+m] = R[k,k]*P[i,j] + T[k,k]*S[i,j]
            end

            UM = SM\[rhs1;rhs2]

            Y[:,k-1] = UM[1:m];Y[:,k]=UM[m+1:end]

            PY[:,k]=P*Y[:,k];PY[:,k-1]=P*Y[:,k-1]
            SY[:,k]=S*Y[:,k];SY[:,k-1]=S*Y[:,k-1]

            k -=2
        end
    end

    if k == 1
        rhs =F[:,1]
        PY[:,2]=P*Y[:,2]
        SY[:,2]=S*Y[:,2]
        for j=2:n
            rhs -= R[1,j]*PY[:,j] + T[1,j]*SY[:,j]
        end

        Y[:,1]=(R[1,1]*P + T[1,1]*S)\rhs
    end
    Y
end

##Solves A*X*transpose(B) + C*X*transpose(D) = E
function lyap(A,B,C,D,E)
    AC=schurfact(Matrix(A),Matrix(C))
    BD=schurfact(Matrix(B),Matrix(D))
    Q1=AC[:left];Q2=BD[:left]
    Z1=AC[:right];Z2=BD[:right]

    Y=lyapuptriang(AC[:S],BD[:S],AC[:T],BD[:T],Q1'*E*Q2)
    Z1*Y*Z2'
end

# function lyap(A,B,F)
#   Λ,V=eig(A)
#   Ω,W=eig(B)
#
#     F=inv(V)*F*W
#
#     Y=lyapdiag!(Λ,Ω,F)
#
#     real(V*Y*inv(W))
# end


text:
################################################################################
#
# Restriction operator linear algebra routines
#
################################################################################

## Restriction operator (explicit matrix)

export restrOpMat

function restrOpMat(n::NTuple{2, Int64}, h::NTuple{2, R}, xint::Array{R, 1}, zint::R, nabc::NTuple{4, Int64})

    # Neighboring indexes for bilinear interp
    idx_x = xint./h[1].+R(1)
    idx_x1 = floor.(idx_x)
    idx_x2 = idx_x1.+R(1)
    idx_z = zint./h[2].+R(1)
    idx_z1 = floor.(idx_z)
    idx_z2 = idx_z1.+R(1)

    j_11 = Int64.([idx_x1...]).+nabc[1].+(Int64(idx_z1)+nabc[3]-1)*(n[1]+nabc[1]+nabc[2])
    j_12 = Int64.([idx_x1...]).+nabc[1].+(Int64(idx_z2)+nabc[3]-1)*(n[1]+nabc[1]+nabc[2])
    j_21 = Int64.([idx_x2...]).+nabc[1].+(Int64(idx_z1)+nabc[3]-1)*(n[1]+nabc[1]+nabc[2])
    j_22 = Int64.([idx_x2...]).+nabc[1].+(Int64(idx_z2)+nabc[3]-1)*(n[1]+nabc[1]+nabc[2])

    # Weights
    v_11 = (1f0.-abs.(idx_x.-idx_x1)).*(1f0.-abs.(idx_z-idx_z1))
    v_12 = (1f0.-abs.(idx_x.-idx_x1)).*(1f0.-abs.(idx_z-idx_z2))
    v_21 = (1f0.-abs.(idx_x.-idx_x2)).*(1f0.-abs.(idx_z-idx_z1))
    v_22 = (1f0.-abs.(idx_x.-idx_x2)).*(1f0.-abs.(idx_z-idx_z2))

    # Assembly
    nx = length(xint)
    i = collect(1:nx)
    next = prod(sz_abc(n, nabc))
    return sparse(i, j_11, v_11, nx, next)+
           sparse(i, j_12, v_12, nx, next)+
           sparse(i, j_21, v_21, nx, next)+
           sparse(i, j_22, v_22, nx, next)

end


text:
using IterativeSolvers
using Test
using Random
using LinearAlgebra
using SparseArrays

@testset "QMR" begin

n = 10
m = 6
Random.seed!(1234567)

@testset "Matrix{$T}" for T in (Float32, Float64, ComplexF32, ComplexF64)
    A = rand(T, n, n) + n * I
    b = rand(T, n)
    tol = √eps(real(T))*10

    x, history = qmr(A, b, log=true)
    @test isa(history, ConvergenceHistory)
    @test history.isconverged
    @test norm(A * x - b) / norm(b) ≤ tol

end

@testset "SparseMatrixCSC{$T, $Ti}" for T in (Float64, ComplexF64), Ti in (Int64, Int32)
    A = sprand(T, n, n, 0.5) + n * I
    b = rand(T, n)
    tol = √eps(real(T))

    x, history = qmr(A, b, log=true)
    @test history.isconverged
    @test norm(A * x - b) / norm(b) ≤ tol
end

@testset "Maximum number of iterations" begin
    x, history = qmr(rand(5, 5), rand(5), log=true, maxiter=2)
    @test history.iters == 2
    @test length(history[:resnorm]) == 2
end

end


text:
using PrettyTables

struct RecurrentRVPTestResults
    results::Vector{RecurrentRVPTestResult}
end

# function RecurrentRVPTestResults(x::MuZeroConfig)
#     getpolicytests()
#     d = policy_testset(x.game_type)
#     tests = []
#     for (idx, kv) in enumerate(d)
#         game, actions = kv
#         push!(tests, RecurrentRVPTest(idx, game, actions))
#     end
#     RecurrentRVPTestResults(tests)
# end

struct RecurrentRVPTestResultsSVG
    tests::RecurrentRVPTestResults
end

using Random

svg(x::RecurrentRVPTestResults) = RecurrentRVPTestResultsSVG(x)
Base.rand(x::RecurrentRVPTestResults) = rand(x.results)
Base.getindex(x::RecurrentRVPTestResults, idx::Integer) = x.results[idx]
Base.getindex(x::RecurrentRVPTestResults, idx::Vector{<:Integer}) = x.results[idx]
Base.length(x::RecurrentRVPTestResults) = length(x.results)
Base.lastindex(x::RecurrentRVPTestResults) = length(x)
Base.iterate(x::RecurrentRVPTestResults, state=1) = state > length(x) ? nothing : (x[state], state+1)

getdf(x::RecurrentRVPTestResults; range=1:length(x), limit=3) = getdf(RecurrentRVPTestResult, x.results[range])

function Base.show(io::IO, m::MIME"text/plain", x::RecurrentRVPTestResults)
    len = length(x.results)
    println()
    pretty_results(
        io,
        getdf(x);
        title="Recurrent RVP Test Results (all $len)"
    )
end

function pretty_results(io::IO, df::DataFrame; title="", pf_col_idx=8)
    pf_hl = Highlighter((data, i, j)->(j==pf_col_idx) && data[i,j] == "Pass", crayon"green");
    pretty_table(
        io,
        df;
        show_row_number=false,
        nosubheader=true,
        title=title,
        highlighters=(pf_hl),
    )
end

function Base.getindex(x::RecurrentRVPTestResults, r::UnitRange{<:Integer})
    println()
    pretty_results(
        stdout,
        getdf(x; range=r),
        title="Policy Test Results ($(r))",
        pf_col_idx=8
    )
end


text:
function Poisson(E, λ)

    (exp.(-λ).*λ.^E)/factorial.(E)

end

function cdfPoisson(x_1, x_2, λ)

    [sum(Poisson.(collect(x_1:i), λ)) for i in x_1:x_2]

end


text:
module AbnormalReturns

using Tables
using LinearAlgebra
using StatsBase
using Reexport
using Statistics
using Dates
using DataFrames
using IntervalSets: ClosedInterval, Ellipsis, (..)
using SparseArrays
@reexport using BusinessDays
@reexport using StatsModels
using ShiftedArrays
using PrettyTables

##############################################################################
##
## Exported methods and types
##
##############################################################################

# types and functions for fast CAR calculations
export MarketData, TimelineTable, car, alpha, beta,
    BasicReg, quick_reg, DataVector,
    bh_return, bhar, MarketCalendar, TimelineColumn

export getindex, values, names, istable, columnaccess, columns,
    getcolumn, columnnames

# From Statistics
export var, std

# From StatsBase
export coef, coefnames, responsename, nobs, dof_residual,
    r2, adjr2, islinear, deviance, rss, predict

export getindex

export dropmissing, select!

export ClosedInterval, ..

export lag, lead

##############################################################################
##
## Load files
##
##############################################################################

include("marketCalendar.jl")
include("dictIndex.jl")
include("timelineData.jl")
include("utils.jl")
include("iterateTimelineTable.jl")
include("fastRegression.jl")
include("calcFunctions.jl")
include("statsModelsModelcols.jl")

end


text:
using Documenter
using QuasiCopula

makedocs(
    sitename = "QuasiCopula",
    format = Documenter.HTML(),
    modules = [QuasiCopula],
    pages = [
        "Home" => "index.md",
        "API" => "man/api.md",
        "AR(1) Covariance" => "man/AR_Examples.md",
        "CS Covariance" => "man/CS_Examples.md",
        "VC Covariance" => "man/VC_Examples.md"
    ]
)

# Documenter can also automatically deploy documentation to gh-pages.
# See "Hosting Documentation" and deploydocs() in the Documenter manual
# for more information.
deploydocs(
    repo   = "github.com/OpenMendel/QuasiCopula.jl.git",
    target = "build"
)


text:
#
#
# Encapulates arbitrary solvers by a debug solver, which is induced by Decoration Options
#
#
export initializeSolver!, doSolverStep!, getSolverResult, stopSolver!
"""
    initializeSolver!(p,o)

Initialize the solver to the optimization [`Problem`](@ref) by initializing
the encapsulated `options` from within the [`RecordOptions`](@ref)` o`.
"""
function initializeSolver!(p::P,o::O) where {P <: Problem, O <: RecordOptions}
    initializeSolver!(p,o.options)
    get(o.recordDictionary, :Start, RecordGroup() )(p, getOptions(o), 0)
    get(o.recordDictionary, :All,   RecordGroup() )(p, getOptions(o), 0)
end
"""
    doSolverStep!(p,o,iter)

Do one iteration step (the `iter`th) for [`Problem`](@ref)` p` by modifying
the values in the [`Options`](@ref)` o.options` and record the result(s).
"""
function doSolverStep!(p::P,o::O, i) where {P <: Problem, O <: RecordOptions}
    doSolverStep!(p,o.options,i)
    get(o.recordDictionary, :Step, RecordGroup() )(p, getOptions(o), i)
    get(o.recordDictionary, :All,   RecordGroup() )(p, getOptions(o), i)
end
"""
    getSolverResult(p,o)

Return the final result after all iterations that is stored within the
(modified during the iterations) [`Options`](@ref)` o`.
"""
function getSolverResult(p::P,o::O) where {P <: Problem, O <: RecordOptions}
    return getSolverResult(p, o.options)
end
"""
    stopSolver!(p,o,i)

determine whether the solver for [`Problem`](@ref) `p` and the
[`RecordOptions`](@ref) `o` should stop at iteration `i`.
If so, do a (final) record to `:All` and `:Stop`.
"""
function stopSolver!(p::P,o::O,i::Int) where {P <: Problem, O <: RecordOptions}
    s = stopSolver!(p,o.options,i)
    if s
        get(o.recordDictionary, :Stop, RecordGroup() )(p, getOptions(o), typemin(Int) )
        get(o.recordDictionary, :All,   RecordGroup() )(p, getOptions(o), typemin(Int ))
    end
    return s
end


text:
@testset "bots" begin

@testset "MCTSBot" begin
    UCT_C = 2.

    init_bot(game, max_simulations, evaluator) = MCTSBot(game, evaluator, UCT_C, max_simulations, 5, true, 42, false, UCT, 0., 0.)

    @testset "can play tic_tac_toe" begin
        game = load_game("tic_tac_toe")
        max_simulations = 100
        evaluator = RandomRolloutEvaluator(20, 42)
        bot0 = init_bot(game, max_simulations, CxxPtr(evaluator))
        bot1 = init_bot(game, max_simulations, CxxPtr(evaluator))
        results = evaluate_bots(new_initial_state(game), [bot0, bot1], 42)
        @test results[1] + results[2] == 0
    end

    @testset "can play single player" begin
        game = load_game("catch")
        max_simulations = 100
        evaluator = RandomRolloutEvaluator(20, 42)
        bot = init_bot(game, max_simulations, CxxPtr(evaluator))
        results = evaluate_bots(new_initial_state(game), [bot], 42)
        @test results[] > 0
    end

    @testset "play three player stochastic games" begin
        game = load_game("pig(players=3,winscore=20,horizon=30)")
        max_simulations = 1000
        evaluator = RandomRolloutEvaluator(20, 42)
        bot0 = init_bot(game, max_simulations, CxxPtr(evaluator))
        bot1 = init_bot(game, max_simulations, CxxPtr(evaluator))
        bot2 = init_bot(game, max_simulations, CxxPtr(evaluator))
        results = evaluate_bots(new_initial_state(game), [bot0, bot1, bot2], 42)
        @test sum(results) == 0
    end

    function get_action_by_str(state, action_str)
        for action in legal_actions(state)
            if action_str == action_to_string(state, current_player(state), action)
                return action
            end
        end
        @error "Illegal action: $action_str"
    end

    function search_tic_tac_toe_state(initial_actions)
        game = load_game("tic_tac_toe")
        state = new_initial_state(game)
        for action_str in split(initial_actions, " ")
            apply_action(state, get_action_by_str(state, action_str))
        end
        evaluator = RandomRolloutEvaluator(20, 42)
        bot = MCTSBot(game, CxxPtr(evaluator), UCT_C, 10000, 10, true, 42, false, UCT, 0., 0.)
        mcts_search(bot, state), state
    end

    @testset "solve draw" begin
        root, state = search_tic_tac_toe_state("x(1,1) o(0,0) x(2,2)")
        @test to_string(state) == "o..\n.x.\n..x"
        @test get_outcome(root)[get_player(root)+1] == 0
        for c in get_children(root)
            @test get_outcome(c)[get_player(c)+1] <= 0
        end
        best = best_child(root)[]
        @test get_outcome(best)[get_player(best)+1] == 0
    end

    @testset "solve loss" begin
        root, state = search_tic_tac_toe_state("x(1,1) o(0,0) x(2,2) o(1,0) x(2,0)")
        @test to_string(state) == "oox\n.x.\n..x"
        @test get_outcome(root)[get_player(root)+1] == -1
        for c in get_children(root)
            @test get_outcome(c)[get_player(c)+1] == -1
        end
    end

    @testset "solve win" begin
        root, state = search_tic_tac_toe_state("x(1,0) o(2,2)")
        @test to_string(state) == ".x.\n...\n..o"
        @test get_outcome(root)[get_player(root)+1] == 1
        best = best_child(root)[]
        @test get_outcome(best)[get_player(best)+1] == 1
        @test action_to_string(state, get_player(best), get_action(best)) == "x(2,0)"
    end
end

end


text:
# Plotting functions

using PyCall, PyPlot
import PyPlot: plot

# Plotting planetary structures

"Label the mass coordinate x-axis"
function label_x_as_mass()
    xlabel(L"Mass enclosed / M$_\oplus$")
end

"Plot the temperature profile of a planet"
function plot_temperature_profile(soln::PlanetStructure{WithTemp}; kws...)
    M = vec(mass(soln))
    T = vec(temperature(soln))
    plot(M/M_earth, T/K; kws...)
    label_x_as_mass()
    ylabel("Temperature / K")
end

"Plot the radius profile of a planet"
function plot_radius_profile(soln::PlanetStructure; kws...)
    M = vec(mass(soln))
    r = vec(radius(soln))
    plot(M/M_earth, r/R_earth; kws...)
    label_x_as_mass()
    ylabel(L"Radius / R$_\oplus$")
end

"Plot the pressure profile of a planet"
function plot_pressure_profile(soln::PlanetStructure; kws...)
    M = vec(mass(soln))
    P = vec(pressure(soln))
    plot(M/M_earth, P/GPa; kws...)
    label_x_as_mass()
    ylabel("Pressure / GPa")
end

"Plot the heat capacity profile of a planet"
function plot_heat_capacity_profile(soln::PlanetStructure{WithTemp},
                                    sys::PlanetSystem{WithTemp}; kws...)
    M = vec(mass(soln))
    P = vec(pressure(soln))
    T = vec(temperature(soln))
    thermal_gradient = sys.structure_equations[3]
    eos = thermal_gradient.eos
    heatcap = thermal_gradient.heatcap

    cₚ = map(heatcap, P, T)
    plot(M/M_earth, cₚ/(J/kg/K); kws...)
    label_x_as_mass()
    ylabel(L"Heat capacity, c$_p$ / J kg$^{-1}$ K$^{-1}$")
end

"Plot the thermal expansivity profile of a planet"
function plot_expansivity_profile(soln::PlanetStructure{WithTemp},
                                  sys::PlanetSystem{WithTemp}; kws...)
    M = vec(mass(soln))
    P = vec(pressure(soln))
    T = vec(temperature(soln))
    thermal_gradient = sys.structure_equations[3]
    eos = thermal_gradient.eos
    heatcap = thermal_gradient.heatcap

    αᵥ = map(M, P, T) do M, P, T
        vs = ValueSet(M, NaN, P, T)
        Ogre.thermexp(eos, vs)
    end

    plot(M/M_earth, αᵥ/(1/K); kws...)
    label_x_as_mass()
    ylabel(L"Thermal expansivity, $\alpha$ / K$^{-1}$")
end

"Plot the density profile of a planet"
function plot_density_profile(soln::PlanetStructure{WithTemp},
                              sys::PlanetSystem{WithTemp}; kws...)
    M = vec(mass(soln))
    P = vec(pressure(soln))
    T = vec(temperature(soln))
    thermal_gradient = sys.structure_equations[3]
    eos = thermal_gradient.eos
    heatcap = thermal_gradient.heatcap

    ρ = map((M, P, T) -> eos(Ogre.ValueSet(M, NaN, P, T)), M, P, T)
    plot(M/M_earth, ρ/(kg/m^3); kws...)
    label_x_as_mass()
    ylabel(L"Density, $\rho$ / kg m$^{-3}$")
end

function plot(soln::PlanetStructure{NoTemp}; kws...)
    fig, axes = subplots(2, 1, sharex=true)

    sca(axes[1])
    plot_pressure_profile(soln)
    setp(axes[1][:get_xticklabels](), visible=false)
    xlabel("")

    sca(axes[2])
    plot_radius_profile(soln)

    for ax in axes
        ax[:yaxis][:get_major_formatter]()[:set_powerlimits]((0, 4))
    end

    tight_layout()
end

function plot(soln::PlanetStructure{WithTemp}; kws...)
    fig, axes = subplots(3, 1, sharex=true)

    sca(axes[1])
    plot_temperature_profile(soln)

    sca(axes[2])
    plot_pressure_profile(soln)

    sca(axes[3])
    plot_radius_profile(soln)

    for ax in axes[1:2]
        setp(ax[:get_xticklabels](), visible=false)
        ax[:set_xlabel]("")
    end

    for ax in axes
        ax[:yaxis][:get_major_formatter]()[:set_powerlimits]((0, 4))
    end

    tight_layout()
end

function plot(soln::PlanetStructure{WithTemp}, sys::PlanetSystem{WithTemp}; kws...)
    fig, axes = subplots(3, 2, sharex=true, figsize=(14,9))

    sca(axes[1])
    plot_temperature_profile(soln)

    sca(axes[2])
    plot_pressure_profile(soln)

    sca(axes[3])
    plot_radius_profile(soln)

    sca(axes[4])
    plot_density_profile(soln, sys)

    sca(axes[5])
    plot_expansivity_profile(soln, sys)

    sca(axes[6])
    plot_heat_capacity_profile(soln, sys)

    for ax in axes[[1,2,4,5]]
        setp(ax[:get_xticklabels](), visible=false)
        ax[:set_xlabel]("")
    end

    for ax in axes
        ax[:yaxis][:get_major_formatter]()[:set_powerlimits]((0, 4))
    end
end

function plot(pb::PhaseBoundary; kwargs...)
    plot(pb.P/Pa, pb.T/K; kwargs...)
end

"Plot the P-T profile of a planet overlaid on the phase boundaries of water"
function phaseplot(soln::PlanetStructure{WithTemp}; kwargs...)
    P = vec(pressure(soln))
    T = vec(temperature(soln))

    figure()

    ax1 = subplot(111)
    xlabel("Pressure / Pa")
    ylabel("Temperature / K")
    xscale("log")
    yscale("log")

    plot_phases()
    plot(P/Pa, T/K, kwargs...)

    xlim(xmin=1e4, xmax=1e12)
    ylim(ymin=200, ymax=10000)

    tight_layout()
end

"Plot the phase boundaries of water in P-T space"
function plot_phases()
    boundaries = WaterData.load_phase_boundaries()["boundaries"]
    iapws = boundaries["iapws"]
    plot(iapws.P/Pa, iapws.T/K, linewidth=1, color="Red")
    for pb in boundaries["dunaeva"]
        plot(pb.P/Pa, pb.T/K, linewidth=1, color="Black")
    end
end


text:
"""
    sortino_ratio(returns; multiplier=1.0, risk_free=0.0)

Calculates the Sortino Ratio, a downside risk-adjusted performance measure. Contrary to the Sharpe Ratio, only deviations below the minimum acceptable returns are included in the calculation of the risk (downside deviation instead of standard deviation).

# Arguments
- `returns`:        Vector of asset returns.
- `multiplier`:     Optional scalar multiplier, i.e. use `12` to annualize monthly returns, and use `252` to annualize daily returns.
- `MAR`:            Optional vector or scalar value denoting the minimum acceptable return(s). Must have same frequency as the provided returns, e.g. daily.

# Sources
- Sortino, F. and Price, L. (1996). Performance Measurement in a Downside Risk Framework. Journal of Investing.
"""
function sortino_ratio(returns; multiplier=1.0, MAR=0.0)
    mean(returns .- MAR) / downside_deviation(returns, MAR; method=:full) * sqrt(multiplier)
end


text:
# print("Running tests...")
@info "Running tests..."
using availability, Test

# @warn "warning !!!"
# @debug "debug !!!"

include("modelTest.jl")
include("serializationTest.jl")

@testset "Availability Calculations" begin
    @test calculate_Avail(supplyTemp, demandTemp) == 90
    @inferred String calculate_Avail(supplyTemp, demandTemp)

    @test_broken calculate_Avail(supplyTemp, demandTemp) == 91
    @test_skip calculate_Avail(supplyTemp, demandTemp) == 91

    # validates the log statements generated by function call
    @test (@test_logs (:info, "calculate availability") calculate_Avail(supplyTemp, demandTemp)) == 90

end

max_stars_repo_name: kawadesumit/availability

max_issues_repo_name: kawadesumit/availability

max_forks_repo_name: kawadesumit/availability


text:
export CUDNNError

struct CUDNNError <: Exception
    code::cudnnStatus_t
end

Base.convert(::Type{cudnnStatus_t}, err::CUDNNError) = err.code

Base.showerror(io::IO, err::CUDNNError) =
    print(io, "CUDNNError: ", name(err), " (code $(reinterpret(Int32, err.code)))")

name(err::CUDNNError) = unsafe_string(cudnnGetErrorString(err))

## API call wrapper

# outlined functionality to avoid GC frame allocation
@noinline function throw_api_error(res)
    throw(CUDNNError(res))
end

function initialize_api()
    CUDA.prepare_cuda_call()
end

macro check(ex)
    quote
        res = @retry_reclaim CUDNN_STATUS_ALLOC_FAILED $(esc(ex))
        if res != CUDNN_STATUS_SUCCESS
            throw_api_error(res)
        end

        return
    end
end


text:
#Load some packages and add workers
using Images
using Distributed

using LinearAlgebra
using Statistics

#Load the package on all workers
addprocs(4)
@everywhere using DPMMSubClusters

#Load image, taken from SINTEL http://sintel.is.tue.mpg.de/
img = load("frame_0001.png")

#Change to channel view
img_channels = channelview(img)
z,x,y = size(img_channels)

#Create input
input_arr = zeros(5,x*y)
for i=1:x
    for j=1:y
        input_arr[:,(i-1)*y+j]=[img_channels[1,i,j],img_channels[2,i,j],img_channels[3,i,j],i,j]
    end
end

#Create HyperParams
#The rgb,xy multiplier allows us to both play with the weight of the xy/rgb
rgb_prior_multiplier = 1
xy_prior_multiplier = 0.1

data_cov = cov(input_arr')
data_cov[4:5,1:3] .= 0
data_cov[1:3,4:5] .= 0

data_cov[1:3,1:3] .*= rgb_prior_multiplier
data_cov[4:5,4:5] .*= xy_prior_multiplier

data_mean = mean(input_arr,dims = 2)[:]

hyper_params = DPMMSubClusters.niw_hyperparams(1.0,
           data_mean,
           8,
           data_cov)

#Run the model
labels,clusters,weights = DPMMSubClusters.fit(input_arr,hyper_params,50000.0,iters = 200, verbose = true)

#Get the cluster color means
color_means = [x.μ[1:3] for x in clusters]

segemnated_image = zeros(3,x,y)
for i=1:x
    for j=1:y
        segemnated_image[:,i,j] = color_means[labels[(i-1)*y+j]]
    end
end
segemnated_image = colorview(RGB,segemnated_image)


text:
function default_scaleselection(X::AbstractMatrix, ::ObsDim.Constant{1})
    1:size(X, 2)
end

function default_scaleselection(X::AbstractMatrix, ::Union{ObsDim.Constant{2},ObsDim.Last})
    1:size(X, 1)
end

function default_scaleselection(x::AbstractVector)
    1:length(x)
end

function default_scaleselection(x::AbstractVector, ::ObsDim.Last)
    1:length(x)
end

function default_scaleselection(x::AbstractVector, ::ObsDim.Constant{M}) where {M}
    collect(1:length(x))
end

function default_scaleselection(D::AbstractDataFrame)
    valid_columns(D)
end

function valid_columns(D::AbstractDataFrame)
    valid_colnames = Symbol[]
    for colname in names(D)
        if (eltype(D[colname]) <: Union{Real,Missing}) & !any(ismissing, D[colname])
            push!(valid_colnames, colname)
        else
            @warn("Skipping \"$colname\" because it either contains missing or is not of type <: Real")
        end
    end
    valid_colnames
end

function valid_columns(D::AbstractDataFrame, colnames)
    valid_colnames = Symbol[]
    for colname in colnames
        if (eltype(D[colname]) <: Union{Real,Missing}) & !(any(ismissing, D[colname]))
            push!(valid_colnames, colname)
        else
            @warn("Skipping \"$colname\" because it either contains missing or is not of type <: Real")
        end
    end
    valid_colnames
end


text:
module DocumenterTools

using DocStringExtensions

"""
    DocumenterTools.generate(path::String = "docs"; name = nothing, format = :html)

Create a documentation stub in `path`, which is usually a sub folder in
the package root. The name of the package is determined automatically,
but can be given with the `name` keyword argument.

`generate` can also be called without any arguments, in which case it simply puts all
the generated files into a `docs` directory in the current working directory. This
way, if you are already in the root directory of your package, you generally only need
to call `generate()` to generate the documentation stub.

`generate` creates the following files in `path`:

```
.gitignore
src/index.md
make.jl
mkdocs.yml
Project.toml
```

# Arguments

**`path`** file path to the documentation directory to be created (default is `"docs"`).

# Keywords Arguments

**`name`** is the name of the package (without `.jl`). If `name` is not given
`generate` tries to detect it automatically.

**`format`** can be either `:html` (default), `:markdown` or `:pdf` corresponding
to the `format` keyword to Documenter's `makedocs` function, see
[Documenter's manual](https://juliadocs.github.io/Documenter.jl/stable/man/other-formats/).

# Examples
```julia-repl
julia> using DocumenterTools

julia> DocumenterTools.generate("path/to/MyPackage/docs")
[ ... output ... ]
```
"""
function generate(path::AbstractString="docs"; name::Union{AbstractString,Nothing}=nothing,
                  format = :html)
    # TODO:
    #   - set up deployment to `gh-pages`
    #   - fetch url and username automatically (e.g from git remote.origin.url)

    if !(format in (:html, :markdown, :pdf))
        throw(ArgumentError("format `$(repr(format))` not supported."))
    end

    path = abspath(path)
    if isdir(path)
        throw(ArgumentError("directory `$(Base.contractuser(path))` already exists."))
    end

    # determine name of package unless it is given
    if name === nothing
        srcdir = normpath(joinpath(path, "..", "src"))
        candidates = String[]
        if isdir(srcdir)
            for file in readdir(srcdir)
                if isfile(joinpath(srcdir, file))
                    modulename = splitext(file)[1]
                    str = read(joinpath(srcdir, file), String)
                    i = findfirst(Regex("module\\s+" * modulename), str)
                    i !== nothing && push!(candidates, modulename)
                end
            end
        end
        if length(candidates) != 1 # || name === nothing
            throw(ArgumentError(string("could not determine name of package located ",
                "in `$(Base.contractuser(normpath(joinpath(path, ".."))))`. ",
                "Please specify the `name` keyword argument to `DocumenterTools.generate`.")))
        end
        name = candidates[1]
        @info("name of package automatically determined to be `$(name)`.")
    end
    @assert name !== nothing

    # deploy the stub
    try
        @info("deploying documentation to `$(Base.contractuser(path))`")
        mkdir(path)

        # create the root doc files
        Generator.savefile(path, ".gitignore") do io
            write(io, Generator.gitignore())
        end
        Generator.savefile(path, "make.jl") do io
            write(io, Generator.make(name; format = format))
        end
        if format === :markdown
            Generator.savefile(path, "mkdocs.yml") do io
                write(io, Generator.mkdocs(name))
            end
        end
        Generator.savefile(path, "Project.toml") do io
            write(io, Generator.project(; format = format))
        end

        # Create the default documentation source files
        Generator.savefile(path, "src/index.md") do io
            write(io, Generator.index(name))
        end
    catch
        rm(path, recursive=true)
        rethrow()
    end
    nothing
end

"""
    DocumenterTools.generate(pkg::Module; dir = "docs", format = :html)

Same as `generate(path::String)` but the `path` and name is determined
automatically from the module.

!!! note
    The package must be in development mode. Make sure you run
    `pkg> develop pkg` from the Pkg REPL, or `Pkg.develop(\"pkg\")`
    before generating docs.

# Examples
```julia-repl
julia> using DocumenterTools

julia> using MyPackage

julia> DocumenterTools.generate(MyPackage)
[ ... output ... ]
```
"""
function generate(pkg::Module; dir = "docs", format = :html)
    package_path = package_devpath(pkg)
    name = String(nameof(pkg))
    return generate(joinpath(package_path, dir); name = name, format = format)
end

"""
$(SIGNATURES)

Returns the path to the top level directory of a devved out package source tree. The package
is identified by its top level module `pkg`.
"""
function package_devpath(pkg::Module)
    pkg == parentmodule(pkg) || throw(ArgumentError("$(pkg) is a submodule. Use the package top-level module."))
    path = pathof(pkg)
    path === nothing && throw(ArgumentError("could not find path to $(pkg)."))
    name = String(nameof(pkg))

    # check that pkg is not originating from a standard installation directory
    # since those are supposed to be immutable.
    for depot in DEPOT_PATH
        sep = Sys.iswindows() ? "\\\\" : "/"
        if startswith(path, joinpath(depot, "packages", name)) &&
            occursin(Regex(name * sep * "\\w{4,5}" * sep * "src" * sep * name * ".jl"), path)
            throw(ArgumentError(string(
                "module $(name) was found in a standard installation directory. ",
                "Please make sure that $(name) is ready for development by running ",
                "`pkg> develop $(name)` from the Pkg REPL, or ",
                "`Pkg.develop(\"$(name)\")` from the Julia REPL, and try again.")))
        end
    end
    # We assume that the path to source file of pkg is ../Package/src/Package.jl, but we
    # return simply the top level directory of the package (i.e. ../Package)
    return normpath(joinpath(path, "..", ".."))
end

include("genkeys.jl")
include("Generator.jl")
include("Themes.jl")

end # module


text:
export Tensor,
       inds,
       store,
       data

# TODO: make into:
# abstract type TensorStorage{El} end <: AbstractVector{El}
abstract type TensorStorage{ElT} end

data(S::TensorStorage) = S.data

Base.@propagate_inbounds Base.getindex(S::TensorStorage,
                                       i::Integer) = getindex(data(S),i)
Base.@propagate_inbounds Base.setindex!(S::TensorStorage,v,
                                        i::Integer) = setindex!(data(S),v,i)

Random.randn!(S::TensorStorage) = randn!(data(S))
Base.fill!(S::TensorStorage,v) = fill!(data(S),v)

Base.convert(::Type{T},D::T) where {T<:TensorStorage} = D

"""
Tensor{StoreT,IndsT}

A plain old tensor (with order independent
interface and no assumption of labels)
"""
struct Tensor{ElT,N,StoreT<:TensorStorage,IndsT} <: AbstractArray{ElT,N}
  store::StoreT
  inds::IndsT
  # The resulting Tensor is a view into the input data
  function Tensor(store::StoreT,inds::IndsT) where {StoreT<:TensorStorage{ElT},IndsT} where {ElT}
    new{ElT,length(IndsT),StoreT,IndsT}(store,inds)
  end
end

store(T::Tensor) = T.store
inds(T::Tensor) = T.inds
ind(T::Tensor,j::Integer) = inds(T)[j]

#
# Tools for working with Dims/Tuples
# TODO: put this in a seperate file
#

# dim and dims are used in the Tensor interface, overload
# base Dims here
dims(ds::Dims) = ds
dense(ds::Dims) = dims(ds)
dense(::Type{DimsT}) where {DimsT<:Dims} = DimsT
dim(ds::Dims) = prod(ds)

Base.length(ds::Type{<:Dims{N}}) where {N} = N

# Used for BlockSparse Tensors
const BlockDims{N} = NTuple{N,Vector{Int}}

Base.length(ds::Type{<:BlockDims{N}}) where {N} = N

# Make the "dense" version of the indices
# For indices with QNs, this means removing the QNs
dense(ds::BlockDims) = dims(ds)
dense(::Type{<:BlockDims{N}}) where {N} = Dims{N}

# This may be a bad idea to overload?
# Type piracy?
Base.strides(ds::Dims) = Base.size_to_strides(1, dims(ds)...)
Base.copy(ds::Dims) = ds

Base.copy(ds::BlockDims) = ds

function dims(ds::BlockDims{N}) where {N}
  return ntuple(i->sum(ds[i]),Val(N))
end
function dim(ds::BlockDims{N}) where {N}
  return prod(dims(ds))
end

# A tuple of the number of blocks in each
# dimension
function nblocks(inds::BlockDims{N}) where {N}
  return ntuple(dim->length(inds[dim]),Val(N))
end

function nblocks(ind::NTuple{N,Int}) where {N}
  return N
end

function blocksize(ind::NTuple{N,Int},i::Int) where {N}
  return ind[i]
end

# Version taking CartestianIndex
function blockdims(inds::BlockDims{N},
                   loc) where {N}
  return ntuple(dim->inds[dim][loc[dim]],Val(N))
end

function whichblock(inds::BlockDims{N},
                    loc::Int) where {N}
  cartesian_loc = CartesianIndices(nblocks(inds))[loc]
  return Tuple(cartesian_loc)
end

# Version taking LinearIndex
function blockdims(inds::BlockDims{N},
                   loc::Int) where {N}
  # TODO: do this without conversion to CartesianIndex?
  # That may involve division and be slow?
  cartesian_loc = CartesianIndices(nblocks(inds))[loc]
  return ntuple(dim->inds[dim][cartesian_loc[dim]],Val(N))
end

function blockdim(inds::BlockDims{N},
                  loc) where {N}
  return prod(blockdims(inds,loc))
end

## TODO: should this be StaticArrays.similar_type?
#Base.promote_rule(::Type{<:Dims},
#                  ::Type{Val{N}}) where {N} = Dims{N}

ValLength(::Type{Dims{N}}) where {N} = Val{N}
ValLength(::Dims{N}) where {N} = Val{N}()

# This is to help with some generic programming in the Tensor
# code (it helps to construct a Tuple(::NTuple{N,Int}) where the
# only known thing for dispatch is a concrete type such
# as Dims{4})
StaticArrays.similar_type(::Type{<:Dims},
                          ::Type{Val{N}}) where {N} = Dims{N}

StaticArrays.similar_type(::Type{<:BlockDims},
                          ::Type{Val{N}}) where {N} = BlockDims{N}

unioninds(is1::Dims{N1},
          is2::Dims{N2}) where {N1,N2} = Dims{N1+N2}((is1...,is2...))

function deleteat(t::NTuple{N},pos::Int) where {N}
  return ntuple(i -> i < pos ? t[i] : t[i+1],Val(N-1))
end

function insertat(t::NTuple{N},
                  val::NTuple{M},
                  pos::Int) where {N,M}
  return ntuple(i -> i < pos ? t[i] :
                ( i > pos+M-1 ? t[i-1] :
                 val[i-pos+1] ), Val(N+M-1))
end

#
# Generic Tensor functions
#

# The size is obtained from the indices
dims(T::Tensor) = dims(inds(T))
dim(T::Tensor) = dim(inds(T))
Base.size(T::Tensor) = dims(T)

Base.copy(T::Tensor) = Tensor(copy(store(T)),copy(inds(T)))

Base.complex(T::Tensor) = Tensor(complex(store(T)),copy(inds(T)))

Random.randn!(T::Tensor) = (randn!(store(T)); return T)

# TODO: For BlockSparse, this needs to include the offsets
function Base.similar(::Type{<:Tensor{ElT,N,StoreT}},dims) where {ElT,N,StoreT}
  return Tensor(similar(StoreT,dim(dims)),dims)
end
# TODO: make sure these are implemented correctly
#Base.similar(T::Type{<:Tensor},::Type{S}) where {S} = Tensor(similar(store(T),S),inds(T))
#Base.similar(T::Type{<:Tensor},::Type{S},dims) where {S} = Tensor(similar(store(T),S),dims)

Base.similar(T::Tensor) = Tensor(similar(store(T)),copy(inds(T)))
Base.similar(T::Tensor,dims) = Tensor(similar(store(T),dim(dims)),dims)
# To handle method ambiguity with AbstractArray
Base.similar(T::Tensor,dims::Dims) = Tensor(similar(store(T),dim(dims)),dims)
Base.similar(T::Tensor,::Type{S}) where {S<:Number} = Tensor(similar(store(T),S),copy(inds(T)))
Base.similar(T::Tensor,::Type{S},dims) where {S<:Number} = Tensor(similar(store(T),S,dim(dims)),dims)
# To handle method ambiguity with AbstractArray
Base.similar(T::Tensor,::Type{S},dims::Dims) where {S<:Number} = Tensor(similar(store(T),S,dim(dims)),dims)

#function Base.convert(::Type{Tensor{<:Number,N,StoreR,Inds}},
#                      T::Tensor{<:Number,N,<:Any,Inds}) where {N,Inds,StoreR}
#  return Tensor(convert(StoreR,store(T)),copy(inds(T)))
#end

function Base.zeros(::Type{<:Tensor{ElT,N,StoreT}},inds) where {ElT,N,StoreT}
  return Tensor(zeros(StoreT,dim(inds)),inds)
end

# This is to fix a method ambiguity with a Base array function
function Base.zeros(::Type{<:Tensor{ElT,N,StoreT}},inds::Dims{N}) where {ElT,N,StoreT}
  return Tensor(zeros(StoreT,dim(inds)),inds)
end

function Base.promote_rule(::Type{<:Tensor{ElT1,N1,StoreT1}},
                           ::Type{<:Tensor{ElT2,N2,StoreT2}}) where {ElT1,ElT2,
                                                                     N1,N2,
                                                                     StoreT1,StoreT2}
  return Tensor{promote_type(ElT1,ElT2),N3,promote_type(StoreT1,StoreT2)} where {N3}
end

function Base.promote_rule(::Type{Tensor{ElT1,N,StoreT1,Inds}},
                           ::Type{Tensor{ElT2,N,StoreT2,Inds}}) where {ElT1,ElT2,N,
                                                                       StoreT1,StoreT2,Inds}
  return Tensor{promote_type(ElT1,ElT2),N,promote_type(StoreT1,StoreT2),Inds}
end

#function Base.promote_rule(::Type{<:Tensor{ElT,<:Any,StoreT,<:Any}},::Type{IndsR}) where {N,ElT,StoreT,IndsR}
#  return Tensor{ElT,length(IndsR),StoreT,IndsR}
#end

# Convert the tensor type to the closest dense
# type
function dense(::Type{<:Tensor{ElT,N,StoreT,IndsT}}) where {ElT,N,StoreT,IndsT}
  return Tensor{ElT,N,dense(StoreT),dense(IndsT)}
end

function StaticArrays.similar_type(::Type{<:Tensor{ElT,<:Any,StoreT,<:Any}},indsR) where {N,ElT,StoreT}
  return Tensor{ElT,length(indsR),StoreT,typeof(indsR)}
end

Base.BroadcastStyle(::Type{T}) where {T<:Tensor} = Broadcast.ArrayStyle{T}()

function Base.similar(bc::Broadcast.Broadcasted{Broadcast.ArrayStyle{T}},
                      ::Type{<:Any}) where {T<:Tensor}
  A = find_tensor(bc)
  return similar(A)
end

# This is used for overloading broadcast
"`A = find_tensor(As)` returns the first Tensor among the arguments."
find_tensor(bc::Base.Broadcast.Broadcasted) = find_tensor(bc.args)
find_tensor(args::Tuple) = find_tensor(find_tensor(args[1]), Base.tail(args))
find_tensor(x) = x
find_tensor(a::Tensor, rest) = a
find_tensor(::Any, rest) = find_tensor(rest)

# TODO: implement some generic fallbacks for necessary parts of the API?
#Base.getindex(A::TensorT, i::Int) where {TensorT<:Tensor} = error("getindex not yet implemented for Tensor type $TensorT")


text:
using Jtac

chain = @chain MetaTac begin
  Conv(128, relu, padding = 1)
  Batchnorm()
#  Conv(128, relu, padding = 1)
#  Batchnorm()
#  Conv(128, relu, padding = 1)
#  Batchnorm()
#  Conv(128, relu, padding = 1)
#  Batchnorm()
#  Conv(128, relu, padding = 1)
#  Batchnorm()
end

vhead = @chain (9, 9, 128) Conv(1, relu, window = 1) Batchnorm() Dense(64, relu) Dense(1)
phead = @chain (9, 9, 128) Conv(2, relu, window = 1) Batchnorm() Dense(81)

#chain = @chain MetaTac Dense(50) Dense(50)
vhead = nothing
phead = nothing

model = NeuralModel(MetaTac, chain, vhead = vhead, phead = phead)
model = model |> to_gpu |> Async

player = MCTSPlayer(model, power = 50, dilution = 0.075)

function train_meta_conv(player)

  mcts_opponents = [MCTSPlayer(power = p, temperature = 0.5) for p in [100, 250, 500, 750]]
  self_opponents = [MCTSPlayer(model, power = p, name = "current$p") for p in [25, 50]]

  loss = Loss(reg = 1e-4)

  i = 0
  for (lr, epochs) in [(2e-2, 20), (1e-3, 40), (2e-4, 60)]
    i += 1
    @time with_contest( train_self!
                , player
                , loss = loss
                , playings = 516
                , batchsize = 1024
                , branching = (before = 0.5, during = 0.0, steps = 1)
                , iterations = 1
                , replays = 2
                , epochs = epochs
                , augment = true
                , opponents = [mcts_opponents; self_opponents]
                , length = 0
                , cache = 2000
                , interval = 10
                , distributed = true
                , optimizer = Knet.Momentum
                , lr = lr
                , gamma = 0.9
                )
    save_model("meta-conv-$i", training_model(player) |> to_cpu)
  end
end

max_stars_repo_name: roSievers/jtac


text:
module YetAnotherHelperFlagKillBox

using ..Ahorn, Maple

@mapdef Trigger "YetAnotherHelper/FlagKillBox" FlagKillBox(x::Integer, y::Integer, width::Integer=Maple.defaultTriggerWidth, height::Integer=Maple.defaultTriggerHeight, flag::String="example_flag")

const placements = Ahorn.PlacementDict(
    "Flag Kill Box (Yet Another Helper)" => Ahorn.EntityPlacement(
        FlagKillBox,
        "rectangle"
    )
)

end

max_stars_repo_name: YetAnotherHelper/YetAnotherHelper

max_issues_repo_name: YetAnotherHelper/YetAnotherHelper

max_forks_repo_name: YetAnotherHelper/YetAnotherHelper


text:
function (+){T<:AbstractFloat, R<:RoundingMode}(a::T, b::T, rounding::R)
    hi,lo = eftAdd(a,b)
    eftRound(hi, lo, rounding)
end

function (-){T<:AbstractFloat, R<:RoundingMode}(a::T, b::T, rounding::R)
    hi,lo = eftSub(a,b)
    eftRound(hi, lo, rounding)
end

function (*){T<:AbstractFloat, R<:RoundingMode}(a::T, b::T, rounding::R)
    hi,lo = eftMul(a,b)
    eftRound(hi, lo, rounding)
end

function (/){T<:AbstractFloat, R<:RoundingMode}(a::T, b::T, rounding::R)
    hi,lo = eftDivApprox(a,b)
    eftRound(hi, lo, rounding)
end

function (sqr){T<:AbstractFloat, R<:RoundingMode}(a::T, rounding::R)
    hi,lo = eftSqr(a)
    eftRound(hi, lo, rounding)
end

function (sqrt){T<:AbstractFloat, R<:RoundingMode}(a::T, rounding::R)
    hi,lo = eftSqrtApprox(a)
    eftRound(hi, lo, rounding)
end

max_stars_repo_name: Jeffrey-Sarnoff/RoundingDirectly

max_issues_repo_name: Jeffrey-Sarnoff/RoundingDirectly

max_forks_repo_name: Jeffrey-Sarnoff/RoundingDirectly


text:
abstract type SaturationGoal end

reached(g::EGraph, goal::Nothing) = false
reached(g::EGraph, goal::SaturationGoal) = false

"""
This goal is reached when the `exprs` list of expressions are in the
same equivalence class.
"""
struct EqualityGoal <: SaturationGoal
    exprs::Vector{Any}
    ids::Vector{Int64}
    function EqualityGoal(exprs, eclasses)
        @assert length(exprs) == length(eclasses) && length(exprs) != 0
        new(exprs, eclasses)
    end
end

function EqualityGoal(g::EGraph, exprs)
    n = length(exprs)
    ids = Vector{Int64}(undef, n)
    for i ∈ 1:n
        ec = addexpr!(g, exprs[i])
        ids[i] = ec.id
    end
    return EqualityGoal(exprs, ids)
end

function reached(g::EGraph, goal::EqualityGoal)
    all(x -> in_same_class(g, goal.ids[1], x), @view goal.ids[2:end])
end

"""
Boolean valued function as an arbitrary saturation goal.
User supplied function must take an [`EGraph`](@ref) as the only parameter.
"""
struct FunctionGoal <: SaturationGoal
    fun::Function
end

function reached(g::EGraph, goal::FunctionGoal)::Bool
    fun(g)
end


text:
using Huggett890, Test

hg = Huggett890;

# function retired_sol_test()
#     @testset "Retired solution" begin
#         sol = init_test_retired_solution();
#         k = 2.3;
#         kp = kprime(sol, k);
#         @test kp isa Float64
#         c = consumption(sol, k);
#         @test c isa Float64
#         @test c > 0.0
#     end
# end

function worker_sol_test()
    @testset "Worker solution" begin
        n = 4;
        sol = init_test_worker_solution(n);
        k = 2.3;
        for eIdx = 1 : n
            kp = kprime(sol, k, eIdx);
            @test kp isa Float64
            c = consumption(sol, k, eIdx);
            @test c isa Float64
            @test c > 0.0

            kpFct = hg.kprime_fct(sol, eIdx);
            kp2 = kpFct(k);
            @test isapprox(kp2, kp)
            cFct = hg.cons_fct(sol, eIdx);
            c2 = cFct(k);
            @test isapprox(c2, c)
        end
    end
end

function solution_test()
    @testset "Solution" begin
        Twork = 4;
        T = 7;
        n = 4;
        k = 3.0;
        sol = init_solution(Twork, T, n);
        @test sol isa Solution

        sol = init_test_solution(Twork, T);
        @test sol isa Solution

        for t = 1 : T
            if hg.isretired(sol, t)
                eIdx = 1;
            else
                eIdx = 3;
            end
            wkSol = hg.get_solution(sol, t);
            kp = hg.kprime(wkSol, k, eIdx);
            kp2 = hg.kprime(sol, t, k, eIdx);
            @test isapprox(kp, kp2)
            c = hg.consumption(wkSol, k, eIdx);
            c2 = consumption(sol, t, k, eIdx);
            @test isapprox(c, c2)
        end
    end
end

function set_solution_test()
    @testset "Set solution" begin
        Twork = 4;
        T = 7;
        n = 4;
        sol = init_solution(Twork, T, n);

        t = 3;
        eIdx = 2;
        wkSol = hg.init_test_worker_solution(n);
        hg.set_solution!(sol, t, wkSol);

        k = 3.4;
        kp = hg.kprime(sol, t, k, eIdx);
        kp2 = hg.kprime(wkSol, k, eIdx);
        @test isapprox(kp, kp2)
    end
end

@testset "Solution" begin
    # retired_sol_test();
    worker_sol_test();
    solution_test();
    set_solution_test();
end

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


text:
using ConfigEnv
using Documenter

DocMeta.setdocmeta!(ConfigEnv, :DocTestSetup, :(using ConfigEnv); recursive=true)

makedocs(;
    modules=[ConfigEnv],
    authors="Andrey Oskin",
    repo="https://github.com/Arkoniak/ConfigEnv.jl/blob/{commit}{path}#{line}",
    sitename="ConfigEnv.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://Arkoniak.github.io/ConfigEnv.jl",
        siteurl="https://github.com/Arkoniak/ConfigEnv.jl",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
        "Usage" => "usage.md",
        "Examples" => "scenarios.md"
    ],
)

deploydocs(;
    repo="github.com/Arkoniak/ConfigEnv.jl",
)


text:
n=100
nr=3
d=randn(n,nr)
g=randn(n,nr)
s=randn(n,nr)
pa=Conv.P_conv(d=d,g=g,s=s, gsize=[n,nr], dsize=[n,nr], ssize=[n,nr]);

function adjtest()
	x=randn(size(F,2))
	y=randn(size(F,1))
	a=LinearAlgebra.dot(y,F*x)
	b=LinearAlgebra.dot(x,adjoint(F)*y)
	c=LinearAlgebra.dot(x, transpose(F)*F*x)
	println("adjoint test: ", a, "\t", b)
	@test isapprox(a,b,rtol=1e-6)
	println("must be positive: ", c)
	@test c>0.0
end

function do_tests(m=3)

	@testset "S" begin
		global pa
		global F=Conv.operator(pa, Conv.S());
		adjtest()
	end

	@testset "G" begin
		global pa
		global F=Conv.operator(pa, Conv.G());
		adjtest()
	end

	@testset "D" begin
		global pa
		global F=Conv.operator(pa, Conv.D());
		adjtest()
	end
end
do_tests()

n=100
nr=3
d=randn(n,nr)
g=randn(n,nr)
s=randn(n)
pa=Conv.P_conv(d=d,g=g,s=s, gsize=[n,nr], dsize=[n,nr], ssize=[n]);
@testset "S" begin
	global pa
	global F=Conv.operator(pa, Conv.S());
	adjtest()
end

#=
@testset "BD" begin
	global p=pa.plsbd
	for attrib in [FBD.S(), FBD.G()]
		global F=FBD.operator(p, attrib);
		adjtest()
	end
end

# adj tests of filter
pa=FBD.P_bandpass(Float64, fmin=0.1, fmax=0.4, nt=101)
global F=FBD.create_operator(pa)
@testset "bandpass" begin
	adjtest()
end
=#


text:
using Test
using Yota
using Yota: gradtape, gradcheck, update_chainrules_primitives!
using Yota: trace, compile, play!
using Yota: make_rrule, YotaRuleConfig
import ChainRulesCore: Tangent, ZeroTangent, rrule_via_ad

# test-only dependencies
using CUDA

include("test_helpers.jl")
include("test_chainrules.jl")
include("test_grad.jl")
include("test_update.jl")
include("test_examples.jl")


text:
# # [Lotka-Volterra Model](@id lotka_volterra_example)
#
#md # [![](https://img.shields.io/badge/show-nbviewer-579ACA.svg)](@__NBVIEWER_ROOT_URL__/examples/lotka-volterra.ipynb)

using Petri
using LabelledArrays
using Plots
using OrdinaryDiffEq

# **Step 1:** Define the states and transitions of the Petri Net
#
# Here we have 2 states, wolves and rabbits, and transitions to
# model predation between the two species in the system

S  = [:rabbits, :wolves]
Δ  = LVector(
       birth=(LVector(rabbits=1), LVector(rabbits=2)),
       predation=(LVector(wolves=1, rabbits=1), LVector(wolves=2)),
       death=(LVector(wolves=1), LVector()),
     )
lotka = Petri.Model(S, Δ)

Graph(lotka)

# **Step 2:** Define the parameters and transition rates
#
# Once a model is defined, we can define out initial parameters `u0`, a time
# span `tspan`, and the transition rates of the interactions `β`

u0 = LVector(wolves=10.0, rabbits=100.0)
tspan = (0.0,100.0)
β = LVector(birth=.3, predation=.015, death=.7);

# **Step 3:** Generate a solver and solve
#
# Finally we can generate a solver and solve the simulation

prob = ODEProblem(lotka, u0, tspan, β)
sol = OrdinaryDiffEq.solve(prob,Tsit5(),reltol=1e-8,abstol=1e-8)

plot(sol)


text:
module CircModDesign
using LinearAlgebra
using PyPlot
# using MPI_Tools
using Interpolations
using Optim
using Elliptic
using Struve
using SpecialFunctions
using FFTW
using ACME
using ElectricalEngineering
using Gtk, DataFrames,PyPlot

μ₀ = 4*π*1e-7

include("Filter_Designer.jl")
include("NoiseModeling.jl")
include("ToroidOptimizer.jl")
include("TransformerDesign.jl")
include("StrayCap.jl")
include("LitzACResist.jl")
include("InductanceCalc.jl")
include("MakeSolenoid.jl")
include("RunAnalysis.jl")
include("SPICE2Matrix.jl")
include("CoupledMagnetics.jl")
include("ThermalModeling.jl")
include("Toroid_Inductance.jl")
include("SkinProxLoss.jl")
include("FieldCalc.jl")
include("Load_LTSpice_Net.jl")
export  PlotImpedanceTransformList,
        BiotSav,
        MakeEllip,
        FieldOnAxis_Circ,
        DesignDriveFilter,
        Par,
        Polar2Cart,
        Z_Cap,
        Z_Ind,
        findResPair,
        Butterworth_2,
        CalcInductance,
        SimpleInduct,
        LitzWire_ACRes,
        MakeEllipSolenoid,
        EvalInduct_Biot,
        EvalField_Centered,
        Noise,
        ThermalNoise,
        CL_DAE_Knight,
        ToroidOptimizer,
        Length2Resist,
        FieldMapPointPath,
        CompareD_Circ_Rect_Toroid,
        GetElementVal,
        UpdateElementVal!,
        UpdateElementESR!,
        SPICE_DF2Matrix_ω,
        PipeFlow,
        OneWireLoss,
        TwoWireLoss,
        MultiWireProximity,
        SeriesWireBundleLosses,
        plot1DMagDiffusion

# Write your package code here.

end


text:
function foo()
    a = 10
    b = 20
    _c1 = (a + b)
    _c2 = (a - b)
    _c3 = (a * b)
    _c4 = (a / b)
    _c5 = -(a)
    d = 2.0
    _e1 = (a + d)
    _e2 = (a - d)
    _e3 = (a * d)
    _e4 = (a / d)
    _f = -3.0
    _g = -(a)
end

function add1(x::Int8, y::Int8)::Int16
    return (x + y)
end

function add2(x::Int16, y::Int16)::Int32
    return (x + y)
end

function add3(x::Int32, y::Int32)::Int64
    return (x + y)
end

function add4(x::Int64, y::Int64)::Int64
    return (x + y)
end

function add5(x::UInt8, y::UInt8)::UInt16
    return (x + y)
end

function add6(x::UInt16, y::UInt16)::UInt32
    return (x + y)
end

function add7(x::UInt32, y::UInt32)::UInt64
    return (x + y)
end

function add8(x::UInt64, y::UInt64)::UInt64
    return (x + y)
end

function add9(x::Int8, y::UInt16)::UInt32
    return (x + y)
end

function sub(x::Int8, y::Int8)::Int8
    return (x - y)
end

function mul(x::Int8, y::Int8)::Int16
    return (x * y)
end

function fadd1(x::Int8, y::Float64)::Float64
    return (x + y)
end

function show()
    rv = fadd1(convert(Int8, 6), 6.0)
    @assert(rv == 12)
    println(join(["OK"], " "))
end

function main()
    foo()
    show()
end

main()


text:
# module escape_dimer_one_hit
# export escape_exit_num_dimer
using DifferentialEquations
using Roots
using Printf
# using MATLAB
# using ProgressMeter

include("leap_frog_definitions.jl")
max_hit_q1 = 1
barrier = 10

condition_max_hits(u, t, integrator) = u[3] > max_hit_q1
affect_stop!(integrator) = terminate!(integrator)
function condition_hits_PSS_q1(u, t, integrator) # Event when event_f(u,t) == 0
    u[1]
end

function affect_update_iterator_q1!(integrator)
    integrator.u[3] = integrator.u[3] + 1
end

callback_max_hits = DiscreteCallback(condition_max_hits, affect_stop!)
callback_hits_PSS_q1 = ContinuousCallback(condition_hits_PSS_q1, nothing, affect_update_iterator_q1!, rootfind = true)

cb = CallbackSet(callback_hits_PSS_q1, callback_max_hits)

function escape_after_one_hit(point, t_end, H)
    Q2 = point[1]; P2 = point[2];
    P1 = P1_find_dimer(Q2, P2, H)
    Q1 = 0

    # if isempty(P1)
    #     P1 = P1_find_dimer_second(Q2, P2, H)
    #     t_end = 1
    # end

    if ~isempty(P1)
        # println(Q1)
        q0, p0 = [zeros(3) for i in 1:2]
        q0[1] = 0
        q0[2] = Q2
        q0[3] = 0
        p0[1] = P1[1]
        p0[2] = P2
        p0[3] = 0
       #constructor for ODE
        prob = HamiltonianProblem{true}(Hamiltonian_Dimer, q0, p0, (0., t_end));
        #solve ode , save_everystep=false is important to prevent sol to include all points, not just event points
        sol=solve(prob, RK4(),maxiters=1e20,callback=cb,save_start=true,save_end=true,save_everystep=false)
        # sol=solve(prob, RK4(),maxiters=1e20,reltol=1e-8,abstol=1e-10,callback=cb,save_start=true,save_end=true,save_everystep=false)
        # sol = solve(prob, Vern9(), maxiters = 1e20, reltol = 1e-13, abstol = 1e-15, callback = cb, save_start = true, save_end = true, save_everystep = false)

        uf = zeros(6)
        uf[1] = sol[1,end]
        uf[2] = sol[2,end]
        uf[3] = sol[3,end]
        uf[4] = sol[4,end]
        uf[5] = sol[5,end]
        uf[6] = sol[6,end]
        if maximum([abs(uf[1]),abs(uf[2]),abs(uf[4]),abs(uf[5])]) > barrier && uf[3] == 1
            q1 = sol[1, 2]
            q2 = sol[2, 2]
            p1 = sol[4, 2]
            p2 = sol[5, 2]
            # q1 = sol[1, end-1]
            # q2 = sol[2, end-1]
            # p1 = sol[4, end-1]
            # p2 = sol[5, end-1]
            q = [q1, q2]
            p = [p1, p2]
            diff = abs(Hamiltonian_Dimer(q, p, 0) - H)

            # prob_sign = HamiltonianProblem{true}(Hamiltonian_Dimer, q0, p0, (0., .001));
            # sol_sign = solve(prob_sign, RK4(), maxiters = 1e20)
            # sign_orientation = sign(sol_sign[1,3])

            if diff < 1e-2
                # return Q1, Q2, P1[1], P2, q1, q2, p1, p2
                return Q1, Q2, P1[1], P2, q1, q2, p1, p2

                # return sol
            end
        end
        # return sol.t, sol[1,:], sol[5,:]
        # else

        # end
        # if uf[3]==1
        #     return uf[1], uf[2], uf[4], uf[5]
        # end

        # if dH<1e-5
        #     if is_it(uf,H,tol_dist)
        #         return sol.u[end][5]
        #     else
        #         return -1
        #     end
        # end1233
    # else
        # return NaN, NaN, NaN, NaN
    end
end

function add_four!(point, radius, one_hit_exit, t_end, energy)
    X = point[2]
    Y = point[4]
    X_plus = X + radius
    Y_plus = Y + radius
    X_minus = X - radius
    Y_minus = Y - radius

    point_2 = [X, Y_plus]

    point_4 = [X_minus, Y]

    point_5 = [X_plus, Y]

    point_7 = [X, Y_minus]

    new_points = [point_2, point_4, point_5, point_7]
    for point in new_points
        new_point = escape_after_one_hit(point, t_end, energy)
        if new_point != nothing
            push!(one_hit_exit, new_point)
        end
    end
end

function add_eight!(point, radius, one_hit_exit, t_end, energy)
    X = point[2]
    Y = point[4]

    X_plus = X + radius
    Y_plus = Y + radius
    X_minus = X - radius
    Y_minus = Y - radius

    point_1 = [X_minus, Y_plus]
    point_2 = [X, Y_plus]
    point_3 = [X_minus, Y_plus]

    point_4 = [X_minus, Y]
    point_5 = [X_plus, Y]

    point_6 = [X_minus, Y_minus]
    point_7 = [X, Y_minus]
    point_8 = [X_minus, Y_minus];

    new_points = [point_1, point_2, point_3, point_4, point_5, point_6, point_7, point_8]
    for point in new_points
        new_point = escape_after_one_hit(point, t_end, energy)
        if new_point != nothing
            push!(one_hit_exit, new_point)
        end
    end
end

function unzipper(array_of_points)
    N = length(array_of_points)
    Q = zeros(N)
    P = zeros(N)
    for (index, point) in enumerate(array_of_points)
        Q[index] = point[1]
        P[index] = point[2]
    end
    return Q, P
end

function unzipper3(array_of_points)
    N=length(array_of_points)
    Q0=zeros(N)
    P0=zeros(N)
    Q1=zeros(N)
    P1=zeros(N)
    for (index,point) in enumerate(array_of_points)
        Q0[index]=point[2]
        P0[index]=point[4]
        Q1[index]=point[6]
        P1[index]=point[8]
    end
    return Q0,P0,Q1,P1
end