text:
#UpperBoundsSymmetric.jl
#symmetric and unimodal assumptions!!!

#generic separation formulation
#E[h(V)] = 0
#H(t) = int_[-t, t] h(1 + s) ds
#sep_fun(a, lam, l, u) optimizes
#   #max a + lam H(t)
#   #s.t. t in [l, u]
#   returns val, tstar
function _vopp_sym_unimodal(S, h, H, sep_fun, N, numCuts, TOL, print_trace)
	#create the symmetric grid
	@assert 1 < S <= 2 "Standardized S must be between 1 and 2 by symmetry"
	ps = collect(range(1, stop=S, length=N + 1))
	ps = vcat(2 .- ps[2:end], ps)
	sort!(ps)  #access by ix + N + 1  so that ix spans from -N:N

	m = Model(solver=GurobiSolver(OutputFlag = false))
	@variable(m, theta)
	@variable(m, lam)
	@variable(m, Q[-N:N] >= 0)
	@constraint(m, sum(Q) == 1)

	#t = 0 case
	@constraint(m, theta + lam * h(1) <= sum(ps[j + N + 1] * Q[j] for j = -N:0))
	@objective(m, Max, theta)

	for iter = 1:numCuts
		status = solve(m)
		feasible = true

		print_trace && println("Iteration: \t", iter)
		for k = 0:N-1
			a = 2 * getvalue(theta)
			a -= 2 * sum(ps[j + N + 1] * getvalue(Q[j]) for j = -N:-(k+1))
			a -= sum(ps[j + N + 1] * getvalue(Q[j]) for j = -k:k)
			rhs = sum(ps[j + N + 1] * getvalue(Q[j]) * (1 - ps[j + N + 1]) for j = -k:k)

			val, tstar = sep_fun(a, getvalue(lam), ps[k + N + 1] - 1,
												   ps[k + 1  + N + 1] - 1)

			if val > rhs + TOL
				feasible = false
				print_trace && println("iter $iter : \t $val \t $rhs")

				@constraint(m, 2tstar * theta + lam * H(tstar) <=
					2tstar * sum(ps[j + N + 1] * Q[j] for j = -N:-(k+1)) +
					sum(ps[j + N + 1] * Q[j] * (1 + tstar - ps[j + N + 1]) for j=-k:k)
					)
			end
		end

		if feasible
			if status != :Optimal
				println("Problem likely infeasible.  Dual Status:\t", status)
				return -1.0
			else #correctly solved
				return 1 / getobjectivevalue(m)
			end
		end
	end
	##Exhausted iteration count
	throw("Max Iteration Count Reached")
	return -1.
end

#exposed function for MAD
function vopp_ub_symmetric_MAD(mu, S, M, D, mode; numCuts=100, N=100, print_trace=false, TOL=1e-6)
	#standardize
	c = mu * (1 - M)
	Sc = vopp.comp_Sc(S, M)
	mode_c = (mode - c) / (mu - c)

	@assert Sc >= 1 "S is too small for given mean.  Problem infeasible."
	if Sc > 2
		#just silently round down for now
		Sc = 2
	end

	#create subfunctions
	h(t) = M * abs(t - 1)/2 - D

	#H(t) = int_[-t, t] h(1 + t)
	function H(t)
		if t < 1e-10
			return 0.
		end
		quadgk(s->h(1+s), -t, 0., t)[1]
	end

	#separator
	#max a + lam H(t)
	#s.t. t in [l, u]
	#returns val, tstar
	function sep_fun(a, lam, l, u)
		val = a * u + lam * H(u)
		tstar = u

		val_l = a * l  + lam * H(l)
		if val_l > val
			val = val_l
			tstar = l
		end

		#compute the middle one if its reasonable.
		if lam != 0
			tmid = 2D / M - a / lam / M
		else
			tmid = -1.
		end

		if l <= tmid <= u
			val_mid = a * tmid + lam * H(tmid)
		else
			val_mid = -Inf
		end

		if val_mid > val
			val = val_mid
			tstar = tmid
		end

		return val, tstar
	end

	#kick it off
	_vopp_sym_unimodal(Sc, h, H, sep_fun, N, numCuts, TOL, print_trace)
end

function vopp_ub_symmetric_CV(mu, S, M, C, mode; TOL=1e-6, N=100, numCuts=500, print_trace=false)
	#standardize
	c = mu * (1 - M)
	Sc = vopp.comp_Sc(S, M)
	mode_c = (mode - c) / (mu - c)

	@assert Sc >= 1 "S is too small for given mean.  Problem infeasible."
	if Sc > 2
		#just silently round down for now
		Sc = 2
	end

	#create subfunctions
	h(t) = M^2 * (t - 1)^2 - C^2

	#H(t) = int_[-t, t] h(1 + t)
	function H(t)
		if t < 1e-10
			return 0.
		end

		quadgk(s->h(1+s), -t, 0, t)[1]
	end

	#separator
	#max a + lam H(t)
	#s.t. t in [l, u]
	#returns val, tstar
	function sep_fun(a, lam, l, u)
		val = a * u + lam * H(u)
		tstar = u

		val_l = a * l + lam * H(l)
		if val_l > val
			val = val_l
			tstar = l
		end

		if lam != 0
			tmid = (2lam * C^2 - a) / (2lam * M^2)
		else
			tmid = -1.
		end

		if tmid > 0
			tmid = sqrt(tmid)
		end
		if l <= tmid <= u
			val_mid = a * tmid + lam * H(tmid)
		else
			val_mid = -Inf
		end

		if val_mid > val
			val = val_mid
			tstar = tmid
		end

		return val, tstar
	end

	#kick it off
	_vopp_sym_unimodal(Sc, h, H, sep_fun, N, numCuts, TOL, print_trace)
end


text:
using FullyThreadedTree
using PrettyTables

max_steps = 12
error_tolerance = 1E-2
R = pi / 12.

# Define a smooth function which has a sharp gradient at some interface
interface(x) = sum(x.^2) - R^2 * (1 + 0.25 * cos(5 * atan(x[end], x[1]) + 1))^2
shape(x) = tanh(32 * pi * interface(x))

function adaptive_refinement(fun::Function, max_steps::Int, error_tolerance; plot_error::Bool = false, print_table::Bool = false, dim::Int = 2)
    tree = Tree(zeros(dim), cell_state = fun)

    collect_data = print_table || plot_error
    if collect_data
        max_error = Vector()
        integral = Vector()
        nr_of_cells = Vector()
        nr_of_active_cells = Vector()
        nr_marked_cells = Vector()
    end

    marked = [tree]
    for steps = 1:max_steps

        refine!(marked, cell_state = fun, recurse = true)
        cells = collect_cells(tree)
        if collect_data
            push!(integral, integrate(tree))
            push!(nr_of_cells, length(cells))
            push!(nr_of_active_cells, length(filter(active, cells)))
            push!(max_error, 0.)
            push!(nr_marked_cells, 0)
        end

        marked = Vector{Tree{dim}}()
        for cell ∈ filter(parent_of_active, cells)
            error_estimate = abs(cell.state - sum([child.state for child ∈ cell.children]) / (1<<dim))
            if collect_data max_error[end] = max(error_estimate, max_error[end]) end
            if error_estimate > error_tolerance
                if collect_data  nr_marked_cells[end] += 1<<dim end
                push!(marked, cell)
            end
        end

        if length(marked) == 0
            max_steps = steps
            break
        end
    end

    if print_table
        formatter = Dict(0 => (v, i) -> typeof(v) == Int ? Int(v) : round(v; digits=5))
        pretty_table(hcat(1 : max_steps, nr_of_active_cells, nr_of_cells, nr_marked_cells, max_error, integral), ["step", "# active_cells", "# cells", "# marked", "error", "integral"], tf = markdown, formatter = formatter)
    end
    if plot_error
        display(plot(nr_of_active_cells, max_error, xaxis = :log, yaxis = :log))
    end

    return tree
end

tree = adaptive_refinement(shape, max_steps, error_tolerance, print_table=true)
display(plot(tree))

ext: jl


text:
using BenchmarkTools

BenchmarkTools.DEFAULT_PARAMETERS.evals = 1;
BenchmarkTools.DEFAULT_PARAMETERS.samples = 1_000;
BenchmarkTools.DEFAULT_PARAMETERS.time_tolerance = 1.0e-8;
BenchmarkTools.DEFAULT_PARAMETERS.overhead = BenchmarkTools.estimate_overhead();

include("thingstotime.jl")

qsys8 = @belapsed sums_toward_half(Rational{Int8},2)
qfir8 = @belapsed sums_toward_half(FastRational{Int8,IsReduced},2)
qfmr8 = @belapsed sums_toward_half(FastRational{Int8,MayReduce},2)

qsys16 = @belapsed sums_toward_half(Rational{Int16},6)
qfir16 = @belapsed sums_toward_half(FastRational{Int16,IsReduced},6)
qfmr16 = @belapsed sums_toward_half(FastRational{Int16,MayReduce},6)

qsys32 = @belapsed sums_toward_half(Rational{Int32},8)
qfir32 = @belapsed sums_toward_half(FastRational{Int32,IsReduced},8)
qfmr32 = @belapsed sums_toward_half(FastRational{Int32,MayReduce},8)

qsys64 = @belapsed sums_toward_half(Rational{Int64},17)
qfir64 = @belapsed sums_toward_half(FastRational{Int64,IsReduced},17)
qfmr64 = @belapsed sums_toward_half(FastRational{Int64,MayReduce},17)

qsys128 = @belapsed sums_toward_half(Rational{Int128},77)
qfir128 = @belapsed sums_toward_half(FastRational{Int128,IsReduced},77)
qfmr128 = @belapsed sums_toward_half(FastRational{Int128,MayReduce},77)

function reltimes(withtype, sys2fir, sys2fmr, fmr2fir)
  print("\n\n$withtype:\t")
  print("sys/fir: $sys2fir\t sys/fmr: $sys2fmr\t fmr/fir: $fmr2fir\n\n")
end

reltimes(Int8, round(qsys8 / qfir8, digits = 2), round(qsys8 / qfmr8, digits = 2), round(qfmr8 / qfir8, digits=3));
reltimes(Int16, round(qsys16 / qfir16, digits = 2), round(qsys16 / qfmr16, digits = 2), round(qfmr16 / qfir16, digits=3));
reltimes(Int32, round(qsys32 / qfir32, digits = 2), round(qsys32 / qfmr32, digits = 2), round(qfmr32 / qfir32, digits=3));
reltimes(Int64, round(qsys64 / qfir64, digits = 2), round(qsys64 / qfmr64, digits = 2), round(qfmr64 / qfir64, digits=3));
reltimes(Int128, round(qsys128 / qfir128, digits = 2), round(qsys128 / qfmr128, digits = 2), round(qfmr128 / qfir128, digits=3));


text:
@testset "Misc/MPolyFactor" begin

  check_fac(esum, a) = check_fac(esum, esum, a)

  function check_fac(low, hi, a)
      f = factor(a)
      g = f.unit
      esum = 0
      for i in f.fac
          g *= i[1]^i[2]
          esum += i[2]
      end
      @test g == a
      @test low <= esum <= hi
  end

  QA, A = PolynomialRing(Hecke.QQ, "A")
  K, a = number_field(A^3 - 2, "a")
  Kxyz, (x, y, z, t) = PolynomialRing(K, ["x", "y", "z", "t"])

  check_fac(0, Kxyz(0))
  check_fac(0, Kxyz(1))
  check_fac(0, Kxyz(2))
  check_fac(6, 2*x^2*y*(1+y)*(x^3-2))
  check_fac(9, 2*x^2*y*(1+y)*(2*x^3-1)*(a*x+a^2*y+1)*(a^2*x+a*y+1)^2)
  check_fac(2, x^6-2*y^3)
  check_fac(2, (3*a*y*z*x^2+z-a)*(7*x+a*z^3*y))
  check_fac(3, (x+y+a)*(x+a*y+1)*(x+y+2*a))
  check_fac(3, (x^2+y^2+a)*(x^3+a*y^3+1)*(x^2+y+2*a))
  check_fac(8, (1+x+a*y+x*y)^3*(1+x+a^2*y+x*y^2)^5)
  check_fac(4, (x-(a+1)^4*y+z)*(a*x+y+z)*(x+a*y+z)*(x+y+a*z))
  check_fac(4, (1+x+y+x^2*(1+z+a*y)+y^2*(z+a*x^2))*
               (a*x+y+z*y)*(1+x+a*y*z+z)*(a+x*y+y+a*z))
  check_fac(2, x^99-2*y^99*z^33)
  check_fac(4, x^66-y^66*z^33)

  function dodd(n::Int)
    count = 0
    for i in 1:2:n
      if n%i == 0
        count += 1
      end
    end
    return count
  end
  for i in 1:10
    check_fac(1 + dodd(i) + Int(i%3 == 0), ((1+x+y+z+t)^i+1)*((1+x+y+z+t)^i+2))
  end

  QA, A = PolynomialRing(Hecke.QQ, "A")
  K,e3 = number_field(A^2+A+1, "e3")
  R,(a,x,y,z,t,u,v,w) = PolynomialRing(K, ["a", "x", "y", "z", "t", "u", "v", "w"])
  check_fac(24,
    3*(x+y)*(x-y)*(x-e3*y)*(x+e3*y)*(x+(-e3-1)*y)*(x+(e3+1)*y)*
    (a^2*y*z^2+e3*x*t*u*v*w)*(a^2*y*u^2+e3*x*z*t*v*w)*(a^2*y*t^2-e3*x*z*u*v*w)*
    (a^2*y*t^2+e3*x*z*u*v*w)*(a^2*y*z^2-e3*x*t*u*v*w)*(a^2*y*u^2-e3*x*z*t*v*w)*
    (a*x*z*t*v+(e3+1)*y*u^2*w^2)*(a*x*t*u*v+(e3+1)*y*z^2*w^2)*
    (a*x*z*t*v+(-e3-1)*y*u^2*w^2)*(a*x*z*u*w+(-e3-1)*y*t^2*v^2)*
    (a*x*z*u*w+(e3+1)*y*t^2*v^2)*(a*x*z*t*w+(e3+1)*y*u^2*v^2)*
    (a*x*z*t*w+(-e3-1)*y*u^2*v^2)*(a*x*z*u*v+(-e3-1)*y*t^2*w^2)*
    (a*x*t*u*w+(e3+1)*y*z^2*v^2)*(a*x*z*u*v+(e3+1)*y*t^2*w^2)*
    (a*x*t*u*v+(-e3-1)*y*z^2*w^2)*(a*x*t*u*w+(-e3-1)*y*z^2*v^2))

end


text:
#Contains all the evaluation stuff
using JLD2
using PyPlot
using FileIO
using Statistics

a = load("Results/11_06_ballbot/1_all.jld2")["result"]
b = load("Results/11_06_ballbot/2_all.jld2")["result"]
c = load("Results/11_06_ballbot/3_all.jld2")["result"]
d = load("Results/11_06_ballbot/4_all.jld2")["result"]
#e = load("Results/Videos/11_06_ballbot/5_all.jld2")["result"]
#d = load("Results/Videos/10_31_swingup/1_all.jld2")["result"]
#e = load("Results/Videos/10_31_swingup/1_all.jld2")["result"]

function f_filt(x)
    l = length(x)
    tau = 0.1
    y = zeros(l)
    y[1] = x[1]

    for i = 2:l
        y[i] = (1-tau)*y[i-1] + tau *x[i]
    end

    return y
end

a_all = 4
iters = 500

RE  = zeros(a_all,iters)
A2C = zeros(a_all,iters)
TNPG= zeros(a_all,iters)
TRPO= zeros(a_all,iters)
PPO = zeros(a_all,iters)
TD3 = zeros(a_all,iters)

RE[1,:]   = a[1,:]
A2C[1,:]  = a[2,:]
TNPG[1,:] = a[3,:]
TRPO[1,:] = a[4,:]
PPO[1,:]  = a[5,:]
TD3[1,:]  = a[6,:]

RE[2,:]   = b[1,:]
A2C[2,:]  = b[2,:]
TNPG[2,:] = b[3,:]
TRPO[2,:] = b[4,:]
PPO[2,:]  = b[5,:]
TD3[2,:]  = b[6,:]

RE[3,:]   = c[1,:]
A2C[3,:]  = c[2,:]
TNPG[3,:] = c[3,:]
TRPO[3,:] = c[4,:]
PPO[3,:]  = c[5,:]
TD3[3,:]  = c[6,:]

RE[4,:]   = d[1,:]
A2C[4,:]  = d[2,:]
TNPG[4,:] = d[3,:]
TRPO[4,:] = d[4,:]
PPO[4,:]  = d[5,:]
TD3[4,:]  = d[6,:]

#=
RE[5,:]   = e[1,:]
A2C[5,:]  = e[2,:]
TNPG[5,:] = e[3,:]
TRPO[5,:] = e[4,:]
PPO[5,:]  = e[5,:]
TD3[5,:]  = e[6,:]
=#

scale  = 1
vis    = 0.5
offset = 0

re_mean = f_filt(vec(mean(RE,dims = 1)) .+ offset)
re_std  = f_filt(scale*vec(std(RE,dims=1)))

a2c_mean = f_filt(vec(mean(A2C,dims = 1)) .+ offset)
a2c_std  = f_filt(scale*vec(std(A2C,dims=1)))

tnpg_mean = f_filt(vec(mean(TNPG,dims = 1)) .+ offset)
tnpg_std  = f_filt(scale*vec(std(TNPG,dims=1)))

trpo_mean = f_filt(vec(mean(TRPO,dims = 1)) .+ offset)
trpo_std  = f_filt(scale*vec(std(TRPO,dims=1)))

ppo_mean = f_filt(vec(mean(PPO,dims = 1)) .+ offset)
ppo_std  = f_filt(scale*vec(std(PPO,dims=1)))

td3_mean = f_filt(vec(mean(TD3,dims = 1)) .+ offset)
td3_std  = f_filt(scale*vec(std(TD3,dims=1)))

figure("ball")
clf()
grid(alpha=0.25)
fill_between(1:iters,re_mean .+ re_std, re_mean .- re_std, alpha=vis)
plot(re_mean, label="REINFORCE")
#fill_between(1:iters,a2c_mean .+ a2c_std, a2c_mean .- a2c_std, alpha=vis)
#plot(a2c_mean, label="A2C")
plot(2,2)
fill_between(1:iters,2,2,2)
#fill_between(1:iters,tnpg_mean .+ tnpg_std, tnpg_mean .- tnpg_std, alpha=vis)
#plot(tnpg_mean, label="TNPG")
plot(2,2)
fill_between(1:iters,2,2,2)
fill_between(1:iters,trpo_mean .+ trpo_std, trpo_mean .- trpo_std, alpha=vis)
plot(trpo_mean, label="TRPO")
fill_between(1:iters,ppo_mean .+ ppo_std, ppo_mean .- ppo_std, alpha=vis)
plot(ppo_mean, label="PPO")
fill_between(1:iters,td3_mean .+ td3_std, td3_mean .- td3_std, alpha=vis)
plot(td3_mean, label="PPO-soft")
xlabel("Iterations", fontsize = 12)
ylabel("Cost", fontsize = 12)
#legend(fontsize = 14)
#ylim([0,15])

RE_min = minimum(RE, dims=2)
RE_min_mean = mean(RE_min)
RE_min_std = std(RE_min)

A2C_min = minimum(A2C, dims=2)
A2C_min_mean = mean(A2C_min)
A2C_min_std = std(A2C_min)

TNPG_min = minimum(TNPG, dims=2)
TNPG_min_mean = mean(TNPG_min)
TNPG_min_std = std(TNPG_min)

TRPO_min = minimum(TRPO, dims=2)
TRPO_min_mean = mean(TRPO_min)
TRPO_min_std = std(TRPO_min)

PPO_min = minimum(PPO, dims=2)
PPO_min_mean = mean(PPO_min)
PPO_min_std = std(PPO_min)

TD3_min = minimum(TD3, dims=2)
TD3_min_mean = mean(TD3_min)
TD3_min_std = std(TD3_min)

max_stars_repo_name: NikonPic/REINFORCEMENT-JULIA

max_issues_repo_name: NikonPic/REINFORCEMENT-JULIA

max_forks_repo_name: NikonPic/REINFORCEMENT-JULIA


text:
# https://adventofcode.com/2020/day/11
input = read("2020/data/day_11.txt", String)

test_input = """
L.LL.LL.LL
LLLLLLL.LL
L.L.L..L..
LLLL.LL.LL
L.LL.LL.LL
L.LLLLL.LL
..L.L.....
LLLLLLLLLL
L.LLLLLL.L
L.LLLLL.LL
"""

function parse_input(input::String)
	lines = map(collect, split(strip(input), "\n"))
    grid = permutedims(hcat(lines...))
	return grid
end

function execute_one_round_part1(grid)
    new = copy(grid)
    nrows, ncols = size(grid)
    for c in 1:ncols
        for r in 1:nrows
            row_range = max(1, r - 1):min(nrows, r + 1)
            col_range = max(1, c - 1):min(ncols, c + 1)
            window = grid[row_range,col_range]

            if grid[r,c] == 'L'  # isempty
                if all(window .∈ ".L")  # all surround empty
                    new[r,c] = '#'
                end
            elseif grid[r,c] == '#' # is occupied
                if sum(window .== '#') >= 5 # 4 or more surround full  (+1 for self)
                    new[r,c] = 'L'
                end
            end
        end
    end
    new
end

function part_1(input)
    grid = parse_input(input)
    old = copy(grid)

    for i in 1:1000000
        new = execute_one_round_part1(old)
        if all(new .== old)
            return sum(new .== '#')
        else
            old = new
        end
    end
end
@assert part_1(test_input) == 37
@info part_1(input)

function get_neighbors(grid, idx)
    all_seats = findall(v -> v == 'L', grid)
    neighbors = []
    nrows, ncols = size(grid)

    # right
    for i in (idx[1] + 1):nrows
        candidate = CartesianIndex(i, idx[2])
        if candidate ∈ all_seats
            push!(neighbors, candidate)
            break
        end
    end

    # bottom right
    c = minimum([nrows - idx[1], ncols - idx[2]])  # ugly
    for i in 1:c
        candidate = CartesianIndex(idx[1] + i, idx[2] + i)
        if candidate ∈ all_seats
            push!(neighbors, candidate)
            break
        end
    end

    # bottom
    for col in (idx[2] + 1):ncols
        candidate = CartesianIndex(idx[1], col)
        if candidate ∈ all_seats
            push!(neighbors, candidate)
            break
        end
    end

    # bottom left
    c = minimum([idx[1] - 1, ncols - idx[2]])  # ugly
    for i in 1:c
        candidate = CartesianIndex(idx[1] - i, idx[2] + i)
        if candidate ∈ all_seats
            push!(neighbors, candidate)
            break
        end
    end

    # left
    for row in (idx[1] - 1):-1:1
        candidate = CartesianIndex(row, idx[2])
        if candidate ∈ all_seats
            push!(neighbors, candidate)
            break
        end
    end

    # top left
    c = minimum([idx[1] - 1, idx[2] - 1])  # ugly
    for i in 1:c
        candidate = CartesianIndex(idx[1] - i, idx[2] - i)
        if candidate ∈ all_seats
            push!(neighbors, candidate)
            break
        end
    end

    # top
    for col in (idx[2] - 1):-1:1
        candidate = CartesianIndex(idx[1], col)
        if candidate ∈ all_seats
            push!(neighbors, candidate)
            break
        end
    end

    # top right
    c = minimum([nrows - idx[1], idx[2] - 1])  # ugly
    for i in 1:c
        candidate = CartesianIndex(idx[1] + i, idx[2] - i)
        if candidate ∈ all_seats
            push!(neighbors, candidate)
            break
        end
    end
    return neighbors
end

function execute_one_round_part2(old, grid)
    new = copy(old)
    nrows, ncols = size(old)
    for c in 1:ncols
        for r in 1:nrows
            neighbors = get_neighbors(grid, CartesianIndex(r, c))
            window = [old[idx] for idx in neighbors]

            if old[r,c] == 'L'  # isempty
                if all(window .∈ ".L")  # all surround empty
                    new[r,c] = '#'
                end
            elseif old[r,c] == '#' # is occupied
                if sum(window .== '#') >= 5 # 4 or more surround full  (+1 for self)
                    new[r,c] = 'L'
                end
            end
        end
    end
    new
end

function part_2(input)
    grid = parse_input(input)
    old = copy(grid)
    for i in 1:1000000
        new = execute_one_round_part2(old, grid)
        if all(new .== old)
            return sum(new .== '#')
        else
            old = new
        end
    end
end
@assert part_2(test_input) == 26
# @info part_2("part 2 $input")

max_stars_repo_name: RensDimmendaal/advent-of-code-julia

max_issues_repo_name: RensDimmendaal/advent-of-code-julia

max_forks_repo_name: RensDimmendaal/advent-of-code-julia


text:
using LinearAlgebra, Distributions
ado_vel_dict = Dict("PEDESTRIAN/1" => MvNormal([0.0, 1.0], Diagonal([0.04, 0.60])));# ado velocity distributions


text:
function f()
    steps = []
    fromLeft = 0
    fromTop = 0
    while !isborder(robot, Nord) || !isborder(robot, West)
        if !isborder(robot, Nord)
            push!(steps, Sud)
            move!(robot, Nord)
            fromTop += 1
        else
            push!(steps, Ost)
            move!(robot, West)
            fromLeft += 1
        end
    end
    i = 0
    while !isborder(robot, Ost)
        if i == fromLeft
            putmarker!(robot)
        end
        move!(robot, Ost)
        i += 1
    end
    i = 0
    while !isborder(robot, Sud)
        if i == fromTop
            putmarker!(robot)
        end
        move!(robot, Sud)
        i += 1
    end
    while !isborder(robot, Nord)
        move!(robot, Nord)
    end
    while !isborder(robot, West)
        move!(robot, West)
    end
    i = 0
    while !isborder(robot, Sud)
        if i == fromTop
            putmarker!(robot)
        end
        move!(robot, Sud)
        i += 1
    end
    i = 0
    while !isborder(robot, Ost)
        if i == fromLeft
            putmarker!(robot)
        end
        move!(robot, Ost)
        i += 1
    end
    while !isborder(robot, West)
        move!(robot, West)
    end
    while !isborder(robot, Nord)
        move!(robot, Nord)
    end
    while length(steps) > 0
        move!(robot, pop!(steps))
    end
end

f()

max_stars_repo_name: TheXlebushek/mirea-progs

max_issues_repo_name: TheXlebushek/mirea-progs

max_forks_repo_name: TheXlebushek/mirea-progs


text:
module KoalaTransforms

# Univariate:
export ChainTransformer
export ToIntTransformer
export UnivariateStandardizer, Standardizer
export OneHotEncoder
export UnivariateBoxCoxTransformer, BoxCoxTransformer
export DataFrameToArrayTransformer, RegressionTargetTransformer
export MakeCategoricalsIntTransformer
export DataFrameToStandardizedArrayTransformer
export IntegerToInt64Transformer
export UnivariateDiscretizer, Discretizer

# extended:
export transform, inverse_transform, fit # from `Koala`

# rexported from Koala:
export FeatureSelector, IdentityTransformer, Machine

# for use in this module:
import Koala: BaseType, params
import Koala: TransformerMachine, Machine
import DataFrames: names, AbstractDataFrame, DataFrame, eltypes
import Distributions
using Statistics

# to be extended:
import Koala: Transformer, transform, fit, inverse_transform
import Base: show

# so above export from Koala works:
import Koala: EmptyTransformer, FeatureSelector
import Koala: IdentityTransformer

# development only:
# using ADBUtilities

# constants:
const N_VALUES_THRESH = 16 # for BoxCoxTransformation

## COMPOSITE TRANSFORMERS

mutable struct ChainTransformer <: Transformer
    chain::Vector{Transformer} # in order applied!
end

# constructor is given individual transforms in order opposite to that
# in which they will be applied:
ChainTransformer(transformers...) = ChainTransformer(reverse(Transformer[transformers...]))

function fit(transformer::ChainTransformer, X::Union{AbstractVector, AbstractDataFrame},
             parallel, verbosity)
    schemes = Array{Any}(undef, length(transformer.chain))

    !isempty(transformer.chain) || error("Empty ChainTransformer.")

    # fit first transformer:
    schemes[1] = fit(transformer.chain[1], X, parallel, verbosity)

    # fit the rest:
    for i in 2:length(transformer.chain)
        X = transform(transformer.chain[i - 1], schemes[i - 1], X)
        schemes[i] = fit(transformer.chain[i], X, parallel, verbosity)
    end

    return schemes
end

function transform(transformer::ChainTransformer, schemes,
                   X::Union{AbstractVector, AbstractDataFrame})
    for i in eachindex(schemes)
        X = transform(transformer.chain[i], schemes[i], X)
    end
    return X
end

function inverse_transform(transformer::ChainTransformer, schemes,
                   X::Union{AbstractVector, AbstractDataFrame})
    for i in reverse(eachindex(schemes))
        X = inverse_transform(transformer.chain[i], schemes[i], X)
    end
    return X
end

function Base.getindex(machine::TransformerMachine{ChainTransformer}, i::Int)
    return TransformerMachine(machine.transformer.chain[end - i + 1],
                              machine.scheme[end - i + 1], :)
end

## FOR RELABELLING BY CONSECUTIVE INTEGERS STARTING AT 1

mutable struct ToIntTransformer <: Transformer
    sorted::Bool
    initial_label::Int # ususally 0 or 1
    map_unseen_to_minus_one::Bool # unseen inputs are transformed to -1
end

ToIntTransformer(; sorted=false, initial_label=1,
                 map_unseen_to_minus_one=false) =
                     ToIntTransformer(sorted, initial_label,
                                      map_unseen_to_minus_one)

struct ToIntScheme{T} <: BaseType
    n_levels::Int
    int_given_T::Dict{T, Int}
    T_given_int::Dict{Int, T}
end

# null scheme constructor:
ToIntScheme(S::Type{T}) where T = ToIntScheme{T}(0, Dict{T, Int}(), Dict{Int, T}())

function fit(transformer::ToIntTransformer, v::AbstractVector{T},
             parallel, verbosity) where T
    int_given_T = Dict{T, Int}()
    T_given_int = Dict{Int, T}()
    vals = collect(Set(v))
    if transformer.sorted
        sort!(vals)
    end
    n_levels = length(vals)
    if n_levels > 2^62 - 1
        error("Cannot encode with integers a vector "*
                         "having more than $(2^62 - 1) values.")
    end
    i = transformer.initial_label
    for c in vals
        int_given_T[c] = i
        T_given_int[i] = c
        i = i + 1
    end
    return ToIntScheme{T}(n_levels, int_given_T, T_given_int)
end

# scalar case:
function transform(transformer::ToIntTransformer, scheme::ToIntScheme{T}, x::T) where T
    ret = 0 # otherwise ret below stays in local scope
    try
        ret = scheme.int_given_T[x]
    catch exception
        if isa(exception, KeyError)
            if transformer.map_unseen_to_minus_one
                ret = -1
            else
                throw(exception)
            end
        end
    end
    return ret
end
inverse_transform(transformer::ToIntTransformer, scheme, y::Int) = scheme.T_given_int[y]

# vector case:
function transform(transformer::ToIntTransformer, scheme::ToIntScheme{T},
                   v::AbstractVector{T}) where T
    return Int[transform(transformer, scheme, x) for x in v]
end
inverse_transform(transformer::ToIntTransformer, scheme::ToIntScheme{T},
                  w::AbstractVector{Int}) where T = T[scheme.T_given_int[y] for y in w]

## Univariate standardization

struct UnivariateStandardizer <: Transformer end

function fit(transformer::UnivariateStandardizer, v::AbstractVector{T},
             parallel, verbosity) where T <: Real
    std(v) > eps(Float64) ||
        @warn "Extremely small standard deviation encountered in standardization."
    return  mean(v), std(v)
end

# for transforming single value:
function transform(transformer::UnivariateStandardizer, scheme, x::Real)
    mu, sigma = scheme
    return (x - mu)/sigma
end

# for transforming vector:
transform(transformer::UnivariateStandardizer, scheme,
          v::AbstractVector{T}) where T <: Real =
              [transform(transformer, scheme, x) for x in v]

# for single values:
function inverse_transform(transformer::UnivariateStandardizer, scheme, y::Real)
    mu, sigma = scheme
    return mu + y*sigma
end

# for vectors:
inverse_transform(transformer::UnivariateStandardizer, scheme,
                  w::AbstractVector{T}) where T <: Real =
                      [inverse_transform(transformer, scheme, y) for y in w]

## Standardization of ordinal features of a DataFrame

mutable struct Standardizer <: Transformer
    features::Vector{Symbol} # features to be standardized; empty means all
end

# lazy keyword constructor:
Standardizer(; features=Symbol[]) = Standardizer(features)

struct StandardizerScheme <: BaseType
    schemes::Matrix{Float64}
    features::Vector{Symbol} # all the feature labels of the data frame fitted
    is_transformed::Vector{Bool}
end

# null scheme:
StandardizerScheme() = StandardizerScheme(zeros(0,0), Symbol[], Bool[])

function fit(transformer::Standardizer, X, parallel, verbosity)

    # determine indices of features to be transformed
    features_to_try = (isempty(transformer.features) ? names(X) : transformer.features)
    is_transformed = Array{Bool}(undef, size(X, 2))
    for j in 1:size(X, 2)
        if names(X)[j] in features_to_try && eltype(X[j]) <: AbstractFloat
            is_transformed[j] = true
        else
            is_transformed[j] = false
        end
    end

    # fit each of those features
    schemes = Array{Float64}(undef, 2, size(X, 2))
    verbosity < 1 || @info "Features standarized: "
    univ_transformer = UnivariateStandardizer()
    for j in 1:size(X, 2)
        if is_transformed[j]
            schemes[:,j] = [fit(univ_transformer, collect(X[j]), true, verbosity - 1)...]
            verbosity < 1 ||
                @info "  :$(names(X)[j])    mu=$(schemes[1,j])  sigma=$(schemes[2,j])"
        else
            schemes[:,j] = Float64[0.0, 0.0]
        end
    end

    return StandardizerScheme(schemes, names(X), is_transformed)

end

function transform(transformer::Standardizer, scheme, X)

    names(X) == scheme.features ||
        error("Attempting to transform data frame with incompatible feature labels.")

    Xnew = X[1:end,:] # make a copy of X, working even for `SubDataFrames`
    univ_transformer = UnivariateStandardizer()
    for j in 1:size(X, 2)
        if scheme.is_transformed[j]
            # extract the (mu, sigma) pair:
            univ_scheme = (scheme.schemes[1,j], scheme.schemes[2,j])
            Xnew[j] = transform(univ_transformer, univ_scheme, collect(X[j]))
        end
    end
    return Xnew

end

## One-hot encoding

"""
    OneHotEncoder(drop_last=false)

Returns a transformer for one-hot encoding the categorical features of
an `AbstractDataFrame` object. Here "categorical" refers to any
feature whose eltype is not a subtype of `AbstractFloat`. All eltypes
must admit a `<` and `string` method.

"""
mutable struct OneHotEncoder <: Transformer
    drop_last::Bool
end

# lazy keyword constructor:
OneHotEncoder(; drop_last::Bool=false) = OneHotEncoder(drop_last)

struct OneHotEncoderScheme <: BaseType
    features::Vector{Symbol}         # feature labels
    spawned_features::Vector{Symbol} # feature labels after one-hot encoding
    values_given_feature::Dict{Symbol,Vector{Any}}
end

function fit(transformer::OneHotEncoder, X::AbstractDataFrame, parallel, verbosity)

    features = names(X)
    values_given_feature = Dict{Symbol, Any}()

    verbosity < 1 || @info "One-hot encoding categorical features."

    for ft in features
        if !(eltype(X[ft]) <: AbstractFloat)
            values_given_feature[ft] = sort!(unique(X[ft]))
            if transformer.drop_last
                values_given_feature[ft] = values_given_feature[ft][1:(end - 1)]
            end
            if verbosity > 1
                n_values = length(keys(values_given_feature[ft]))
                @info "Spawned $n_values columns to one-hot encode $ft."
            end
        end
    end

    spawned_features = Symbol[]

    for ft in features
        if !(eltype(X[ft]) <: AbstractFloat)
            for value in values_given_feature[ft]
                subft = Symbol(string(ft,"__",value))

                # in the (rare) case subft is not a new feature label:
                while subft in features
                    subft = Symbol(string(subft,"_"))
                end

                push!(spawned_features, subft)
            end
        else
            push!(spawned_features, ft)
        end
    end

    return OneHotEncoderScheme(features, spawned_features, values_given_feature)

end

function transform(transformer::OneHotEncoder, scheme, X::AbstractDataFrame)

    Set(names(X)) == Set(scheme.features) ||
        error("Attempting to transform DataFrame with incompatible feature labels.")

    # todo: check matching eltypes

    Xout = DataFrame()
    for ft in scheme.features
        if !(eltype(X[ft]) <: AbstractFloat)
            for value in scheme.values_given_feature[ft]
                subft = Symbol(string(ft,"__",value))

                # in case subft is not a new feature name:
                while subft in scheme.features
                    subft = Symbol(string(subft,"_"))
                end

                subft_col = map(X[ft]) do x
                    x == value ? 1.0 : 0.0
                end
                Xout[subft] = convert(Array{Float64}, subft_col)
            end
        else
            Xout[ft] = X[ft]
        end
    end
    return Xout

end

## Univariate Box-Cox transformations

function standardize(v)
    map(v) do x
        (x - mean(v))/std(v)
    end
end

function midpoints(v::AbstractVector{T}) where T <: Real
    return [0.5*(v[i] + v[i + 1]) for i in 1:(length(v) -1)]
end

function normality(v)

    n  = length(v)
    v = standardize(convert(Vector{Float64}, v))

    # sort and replace with midpoints
    v = midpoints(sort!(v))

    # find the (approximate) expected value of the size (n-1)-ordered statistics for
    # standard normal:
    d = Distributions.Normal(0,1)
    w= map(collect(1:(n-1))/n) do x
        quantile(d, x)
    end

    return cor(v, w)

end

function boxcox(lambda, c, x::Real)
    c + x >= 0 || throw(DomainError)
    if lambda == 0.0
        c + x > 0 || throw(DomainError)
        return log(c + x)
    end
    return ((c + x)^lambda - 1)/lambda
end

boxcox(lambda, c, v::AbstractVector{T}) where T <: Real =
    [boxcox(lambda, c, x) for x in v]

"""
## `struct UnivariateBoxCoxTransformer`

A type for encoding a Box-Cox transformation of a single variable
taking non-negative values, with a possible preliminary shift. Such a
transformation is of the form

    x -> ((x + c)^λ - 1)/λ for λ not 0
    x -> log(x + c) for λ = 0

###  Usage

    `trf = UnivariateBoxCoxTransformer(; n=171, shift=false)`

Returns transformer that on fitting to data (see below) will try `n`
different values of the Box-Cox exponent λ (between `-0.4` and `3`) to
find an optimal value. If `shift=true` and zero values are encountered
in the data then the transformation sought includes a preliminary
positive shift by `0.2` times the data mean. If there are no zero
values, then `s.c=0`.

See also `BoxCoxScheme` a transformer for selected ordinals in a DataFrame.

"""
mutable struct UnivariateBoxCoxTransformer <: Transformer
    n::Int      # nbr values tried in optimizing exponent lambda
    shift::Bool # whether to shift data away from zero
end

# lazy keyword constructor:
UnivariateBoxCoxTransformer(; n=171, shift=false) = UnivariateBoxCoxTransformer(n, shift)

function fit(transformer::UnivariateBoxCoxTransformer, v::AbstractVector{T},
    parallel, verbosity) where T <: Real

    m = minimum(v)
    m >= 0 || error("Cannot perform a Box-Cox transformation on negative data.")

    c = 0.0 # default
    if transformer.shift
        if m == 0
            c = 0.2*mean(v)
        end
    else
        m != 0 || error("Zero value encountered in data being Box-Cox transformed.\n"*
                        "Consider calling `fit!` with `shift=true`.")
    end

    lambdas = range(-0.4, stop=3, length=transformer.n)
    scores = Float64[normality(boxcox(l, c, v)) for l in lambdas]
    lambda = lambdas[argmax(scores)]

    return  lambda, c

end

# for X scalar or vector:
transform(transformer::UnivariateBoxCoxTransformer, scheme, X) =
    boxcox(scheme..., X)

# scalar case:
function inverse_transform(transformer::UnivariateBoxCoxTransformer,
                           scheme, x::Real)
    lambda, c = scheme
    if lambda == 0
        return exp(x) - c
    else
        return (lambda*x + 1)^(1/lambda) - c
    end
end

# vector case:
function inverse_transform(transformer::UnivariateBoxCoxTransformer,
                           scheme, w::AbstractVector{T}) where T <: Real
    return [inverse_transform(transformer, scheme, y) for y in w]
end

## `BoxCoxScheme`

"""
## `struct BoxCoxTransformer`

Transformer for Box-Cox transformations on each numerical feature of a
`DataFrame` object. Here "numerical" means of eltype `T <:
AbstractFloat`.

### Method calls

To calculate the compute Box-Cox transformations of a DataFrame `X`
and transform a new DataFrame `Y` according to the same
transformations:

    julia> transformer = BoxCoxTransformer()
    julia> transformerM = TransformerMachine(transformer, X)
    julia> transform(transformerM, Y)

### Transformer parameters

Calls to the first method above may be issued with the following
keyword arguments, with defaluts as indicated:

- `shift=true`: allow data shift in case of fields taking zero values
(otherwise no transformation will be applied).

- `n=171`: number of values of exponent `lambda` to try during optimization.

## See also

`UnivariateBoxCoxTransformer`: The single variable version of the same transformer.

"""
mutable struct BoxCoxTransformer <: Transformer
    n::Int                     # number of values considered in exponent optimizations
    shift::Bool                # whether or not to shift features taking zero as value
    features::Vector{Symbol}   # features to attempt fitting a
                               # transformation to (empty means all)
end

# lazy keyword constructor:
BoxCoxTransformer(; n=171, shift = false, features=Symbol[]) =
    BoxCoxTransformer(n, shift, features)

struct BoxCoxTransformerScheme <: BaseType
    schemes::Matrix{Float64} # each col is a [lambda, c]' pair; one col per feature
    features::Vector{Symbol} # all features in the dataframe that was fit
    feature_is_transformed::Vector{Bool} # keep track of which features are transformed
end

# null scheme:
BoxCoxTransformerScheme() = BoxCoxTransformerScheme(zeros(0,0),Symbol[],Bool[])

function fit(transformer::BoxCoxTransformer, X, parallel, verbosity)

    verbosity < 1 || @info "Computing Box-Cox transformations "*
                          "on numerical features."
    # determine indices of features to be transformed
    features_to_try = (isempty(transformer.features) ? names(X) : transformer.features)
    feature_is_transformed = Array{Bool}(undef, size(X, 2))
    for j in 1:size(X, 2)
        if names(X)[j] in features_to_try && eltype(X[j]) <:
            AbstractFloat && minimum(X[j]) >= 0
            feature_is_transformed[j] = true
        else
            feature_is_transformed[j] = false
        end
    end

    # fit each of those features with best Box Cox transformation
    schemes = Array{Float64}(undef, 2, size(X,2))
    univ_transformer = UnivariateBoxCoxTransformer(shift=transformer.shift,
                                               n=transformer.n)
    verbosity < 2 ||
        @info "Box-Cox transformations: "
    for j in 1:size(X,2)
        if feature_is_transformed[j]
            if minimum(X[j]) == 0 && !transformer.shift
                verbosity < 2 ||
                    @info "  :$(names(X)[j])    "*
                            "(*not* transformed, contains zero values)"
                feature_is_transformed[j] = false
                schemes[:,j] = [0.0, 0.0]
            else
                n_values = length(unique(X[j]))
                if n_values < N_VALUES_THRESH
                    verbosity < 2 ||
                        @info "  :$(names(X)[j])    "*
                                "(*not* transformed, less than $N_VALUES_THRESH values)"
                    feature_is_transformed[j] = false
                    schemes[:,j] = [0.0, 0.0]
                else
                    lambda, c = fit(univ_transformer, collect(X[j]), true, verbosity-1)
                    if lambda in [-0.4, 3]
                        verbosity < 2 ||
                            @info "  :$(names(X)[j])    "*
                                    "(*not* transformed, lambda too extreme)"
                        feature_is_transformed[j] = false
                        schemes[:,j] = [0.0, 0.0]
                    elseif lambda == 1.0
                        verbosity < 2 ||
                            @info "  :$(names(X)[j])    "*
                                    "(*not* transformed, not skewed)"
                        feature_is_transformed[j] = false
                        schemes[:,j] = [0.0, 0.0]
                    else
                        schemes[:,j] = [lambda, c]
                        verbosity <1 ||
                            @info "  :$(names(X)[j])    lambda=$lambda  "*
                                    "shift=$c"
                    end
                end
            end
        else
            schemes[:,j] = [0.0, 0.0]
        end
    end

    if !transformer.shift && verbosity < 1
        @info "To transform non-negative features with zero values use shift=true."
    end

    return BoxCoxTransformerScheme(schemes, names(X), feature_is_transformed)

end

function transform(transformer::BoxCoxTransformer, scheme, X::AbstractDataFrame)

    names(X) == scheme.features ||
        error("Attempting to transform a data frame with  "*
              "incompatible feature labels.")

    univ_transformer = UnivariateBoxCoxTransformer(shift=transformer.shift,
                                               n=transformer.n)

    Xnew = X[1:end,:] # make a copy of X
    for j in 1:size(X, 2)
        if scheme.feature_is_transformed[j]
            try
                # extract the (lambda, c) pair:
                univ_scheme = (scheme.schemes[1,j], scheme.schemes[2,j])

                Xnew[j] = transform(univ_transformer, univ_scheme, collect(X[j]))
            catch DomainError
                @warn "Data outside of the domain of the fitted Box-Cox"*
                      " transformation scheme encountered in feature "*
                      "$(names(df)[j]). Transformed to zero."
            end
        end
    end

    return Xnew

end

## General purpose transformer for supervised learning algorithms
## needing floating point matrix input.

"""
    DataFrameToArrayTransformer(boxcox=false, shift=true,
    standardize=false, drop_last=false)

Use this transformer to prepare dataframe inputs for use with
supervised learning algorithms requiring `Matrix{T}` inputs, for
`T<:AbstractFloat`. Here `T` is forced to be `Float64`. Includes one-hot
encoding of categoricals (any feature not of eltype T<:AbstractFloat) and an
option for Box-Cox transformations. Note that all eltypes must admit
`<` and `string` methods.

### Keyword arguments:

- `boxcox`:  whether to apply Box-Cox transformations to the columns

- `shift`: do we shift away from zero in Box-Cox transformations?

- `standardize`: whether to standardize the columns

- `drop_last`: do we drop the last slot in the one-hot-encoding?

"""
mutable struct DataFrameToArrayTransformer <: Transformer
    boxcox::Bool
    shift::Bool
    standardize::Bool
    drop_last::Bool
end

# lazy keyword constructors:
DataFrameToArrayTransformer(; boxcox=false, shift=true, standardize=false,
                            drop_last=false) =
    DataFrameToArrayTransformer(boxcox, shift, standardize, drop_last)

struct Scheme_X <: BaseType
    boxcox::BoxCoxTransformerScheme
    stand::StandardizerScheme
    hot::OneHotEncoderScheme
    features::Vector{Symbol}
    spawned_features::Vector{Symbol} # ie after one-hot encoding
end

function fit(transformer::DataFrameToArrayTransformer, X::AbstractDataFrame, parallel, verbosity)

    features = names(X)

    # fit Box-Cox transformation to numerical features:
    if transformer.boxcox
        verbosity < 1 || @info "Determining Box-Cox transformation parameters."
        boxcox_transformer = BoxCoxTransformer(shift=transformer.shift)
        boxcox = fit(boxcox_transformer, X, true, verbosity - 1)
        X = transform(boxcox_transformer, boxcox, X)
    else
        boxcox = BoxCoxTransformerScheme() # null scheme
    end

    if transformer.standardize
        verbosity < 1 || @info "Determining standardization parameters."
        standardizer = Standardizer()
        stand = fit(standardizer, X, true, verbosity - 1)
        X = transform(standardizer, stand, X)
    else
        stand = StandardizerScheme() # null scheme
    end

    verbosity < 1 || @info "Determining one-hot encodings for data frame categoricals."
    hot_transformer = OneHotEncoder(drop_last=transformer.drop_last)
    hot =  fit(hot_transformer, X, true, verbosity - 1)
    spawned_features = hot.spawned_features

    return Scheme_X(boxcox, stand, hot, features, spawned_features)

end

function transform(transformer::DataFrameToArrayTransformer, scheme_X, X::AbstractDataFrame)
    issubset(Set(scheme_X.features), Set(names(X))) ||
        error("DataFrame feature incompatibility encountered.")
    X = X[scheme_X.features]

    if transformer.boxcox
        boxcox_transformer = BoxCoxTransformer(shift=transformer.shift)
        X = transform(boxcox_transformer, scheme_X.boxcox, X)
    end

    if transformer.standardize
        X = transform(Standardizer(), scheme_X.stand, X)
    end

    hot_transformer = OneHotEncoder(drop_last=transformer.drop_last)
    X = transform(hot_transformer, scheme_X.hot, X)
    return convert(Array{Float64}, X)
end

"""
## `RegressionTargetTransformer(boxcox=false, shift=true, standardize=true)::Transformer`

A general purpose transformer for target variables in regression
problems. Standardizes by default, and includes Box-Cox
transformations as an option.

### Keyword arguments:

- `boxcox`:  whether to apply Box-Cox transformations to the input patterns

- `shift`: do we shift away from zero in Box-Cox transformations?

- `standardize`: do we standardize (after any Box-Cox transformations)?

"""
mutable struct RegressionTargetTransformer <: Transformer
    boxcox::Bool # do we apply Box-Cox transforms to target (before any standarization)?
    shift::Bool # do we shift away from zero in Box-Cox transformations?
    standardize::Bool # do we standardize targets?
end

# lazy keyword constructors:
RegressionTargetTransformer(; boxcox=false, shift=true, standardize=true) =
    RegressionTargetTransformer(boxcox, shift, standardize)

struct Scheme_y <: BaseType
    boxcox::Tuple{Float64,Float64}
    standard::Tuple{Float64,Float64}
end

function fit(transformer::RegressionTargetTransformer, y, parallel, verbosity)

    if transformer.boxcox
        verbosity < 1 || @info "Computing Box-Cox transformations for target."
        boxcox_transformer = UnivariateBoxCoxTransformer(shift=transformer.shift)
        boxcox = fit(boxcox_transformer, y, true, verbosity - 1)
        y = transform(boxcox_transformer, boxcox, y)
    else
        boxcox = (0.0, 0.0) # null scheme
    end
    if transformer.standardize
        verbosity < 1 || @info "Computing target standardization."
        standard_transformer = UnivariateStandardizer()
        standard = fit(standard_transformer, y, true, verbosity - 1)
    else
        standard = (0.0, 1.0) # null scheme
    end
    return Scheme_y(boxcox, standard)
end

function transform(transformer::RegressionTargetTransformer, scheme_y, y)
    if transformer.boxcox
        boxcox_transformer = UnivariateBoxCoxTransformer(shift=transformer.shift)
        y = transform(boxcox_transformer, scheme_y.boxcox, y)
    end
    if transformer.standardize
        standard_transformer = UnivariateStandardizer()
        y = transform(standard_transformer, scheme_y.standard, y)
    end
    return y
end

function inverse_transform(transformer::RegressionTargetTransformer, scheme_y, y)
    if transformer.standardize
        standard_transformer = UnivariateStandardizer()
        y = inverse_transform(standard_transformer, scheme_y.standard, y)
    end
    if transformer.boxcox
        boxcox_transformer = UnivariateBoxCoxTransformer(shift=transformer.shift)
        y = inverse_transform(boxcox_transformer, scheme_y.boxcox, y)
    end
    return y
end

## TRANSFORMER TO CONVERT CATEGORICALS TO INTEGER-VALUED COLUMNS

mutable struct MakeCategoricalsIntTransformer <: Transformer
    initial_label::Int
    sorted::Bool
end

MakeCategoricalsIntTransformer(; initial_label=1, sorted=false) = MakeCategoricalsIntTransformer(initial_label, sorted)

struct MakeCategoricalsIntScheme <: BaseType
    categorical_features::Vector{Symbol}
    schemes::Vector{ToIntScheme}
    to_int_transformer::ToIntTransformer
end

function fit(transformer::MakeCategoricalsIntTransformer, X::AbstractDataFrame, parallel, verbosity)
    categorical_features = Symbol[]
    schemes = ToIntScheme[]
    to_int_transformer = ToIntTransformer(sorted=transformer.sorted,
                                          initial_label=transformer.initial_label)
    types = eltypes(X)
    features = names(X)
    for j in eachindex(types)
        if !(types[j] <: AbstractFloat)
            push!(categorical_features, features[j])
            push!(schemes, fit(to_int_transformer, X[j], parallel, verbosity))
        end
    end
    verbosity < 1 || @info "Input features treated as categorical: $categorical_features"
    return MakeCategoricalsIntScheme(categorical_features, schemes, to_int_transformer)
end

function transform(transformer::MakeCategoricalsIntTransformer, scheme_X, X::AbstractDataFrame)
    Xt = X[1:end,:] # make a copy of X, working even for `SubDataFrame`s
    for j in eachindex(scheme_X.categorical_features)
        ftr = scheme_X.categorical_features[j]
        Xt[ftr] = transform(scheme_X.to_int_transformer, scheme_X.schemes[j], X[ftr])
    end
    return Xt
end

## DISCRETIZATION OF CONTINUOUS VARIABLES

"""
    UnivariateDiscretizer(n_classes=512)

Returns a `Koala.Transformer` for discretizing vectors of `Real`
eltype, where `n_classes` describes the resolution of the
discretization. Transformed vectors are of eltype `Int46`. The
transformation is chosen so that the vector on which the transformer
is fit has, in transformed form, an approximately uniform distribution
of values.

### Example

    using Koala
    using KoalaTransforms
    t = UnivariateDiscretizer(n_classes=10)
    v = randn(1000)
    tM = Machine(t, v)   # fit the transformer on `v`
    w = transform(tM, v) # transform `v` according to `tM`
    StatsBase.countmap(w)

"""
mutable struct UnivariateDiscretizer <: Transformer
    n_classes::Int
end

# lazy keyword constructor:
UnivariateDiscretizer(; n_classes=512) = UnivariateDiscretizer(n_classes)

struct UnivariateDiscretizerScheme <: BaseType
    odd_quantiles::Vector{Float64}
    even_quantiles::Vector{Float64}
end

function fit(transformer::UnivariateDiscretizer, v, parallel, verbosity)
    n_classes = transformer.n_classes
    quantiles = quantile(v, Array(range(0, stop=1, length=2*n_classes+1)))
    clipped_quantiles = quantiles[2:2*n_classes] # drop 0% and 100% quantiles

    # odd_quantiles for transforming, even_quantiles used for inverse_transforming:
    odd_quantiles = clipped_quantiles[2:2:(2*n_classes-2)]
    even_quantiles = clipped_quantiles[1:2:(2*n_classes-1)]

    return UnivariateDiscretizerScheme(odd_quantiles, even_quantiles)
end

# transforming scalars:
function transform(transformer::UnivariateDiscretizer, scheme, r::Real)
    k = 1
    for level in scheme.odd_quantiles
        if r > level
            k = k + 1
        end
    end
    return k
end

# transforming vectors:
function transform(transformer::UnivariateDiscretizer, scheme,
                   v::AbstractVector{T}) where T<:Real
    return [transform(transformer, scheme, r) for r in v]
end

# scalars:
function inverse_transform(transformer::UnivariateDiscretizer, scheme, k::Int)
    n_classes = length(scheme.even_quantiles)
    if k < 1
        return scheme.even_quantiles[1]
    elseif k > n_classes
        return scheme.even_quantiles[n_classes]
    end
    return scheme.even_quantiles[k]
end

# vectors:
function inverse_transform(transformer::UnivariateDiscretizer, scheme,
                           w::AbstractVector{T}) where T<:Integer
    return [inverse_transform(transformer, scheme, k) for k in w]
end

## CONVERTING ANY INTEGER TYPE TO INT16

mutable struct IntegerToInt64Transformer <: Transformer
end

fit(transformer::IntegerToInt64Transformer, X, parallel, verbosity) = nothing

transform(transformer::IntegerToInt64Transformer, scheme, X) = Int64[X...]

## DISCRETIZING ALL COLUMNS OF A DATAFRAME

struct NominalOrdinalIntArray <: BaseType
    A::Matrix{Int}
    features::Vector{Symbol}
    is_ordinal::Vector{Bool}
end

function show(stream::IO, ::MIME"text/plain", X::NominalOrdinalIntArray)
    features_plus = Array{String}(undef, length(X.features))
    for j in eachindex(X.features)
        kind = X.is_ordinal[j] ? "ordinal" : "nominal"
        features_plus[j] = string(X.features[j], " ($kind)  ")
    end
    for str in features_plus
        print(stream, str)
    end
    println(stream)
    show(stream, MIME("text/plain"), X.A)
end

"""
    Discretizer(n_classes=512)

Returns a `Koala.Transformer` for discretizing the columns of a
`DataFrame`. Transformed objects are of type `NominalOrdinalIntArray`,
a lightweight wrapper for an `Int64` matrix whose fields `A`,
`features` and `is_ordinal` store, respectively, the integer matrix, a
vector of feature names, and a vector recording which features are
ordinal.

Columns are discretized according to the following conventions: If the
column eltype is a subtype of `AbstractFloat`, the column is
discretized using a `UnivariateDiscretizer` transformer. All remaining
eltypes are treated as nominal, with elements relabeled with integers,
starting with 1; nominal values not encountered during fitting the
transformer will throw a `KeyError` if there is an attempt to
transform them according to the fitted scheme.

"""
mutable struct Discretizer <: Transformer
    n_classes::Int
    features::Vector{Symbol} # features to be used; empty means all
end

# lazy keyword constructor:
Discretizer(; n_classes=512, features=Symbol[]) = Discretizer(n_classes, features)

struct DiscretizerScheme
    features::Vector{Symbol} # features actually used
    transformer_machines::Vector{TransformerMachine}
    is_ordinal::Vector{Bool}
end

function fit(transformer::Discretizer, X::AbstractDataFrame, parallel, verbosity)

    to_int = ToIntTransformer(sorted=true)
    discrete = UnivariateDiscretizer()

    if isempty(transformer.features)
        features = names(X)
    else
        features = transformer.features
    end

    transformer_machines = Array{TransformerMachine}(undef, length(features))
    is_ordinal = Array{Int}(undef, length(features))

    j = 1
    for ftr in features
        if eltype(X[ftr]) <: AbstractFloat
            is_ordinal[j] = true
            transformer_machines[j] = Machine(discrete, X[ftr];
                parallel=false, verbosity=verbosity - 1)
        else
            is_ordinal[j] = false
            transformer_machines[j] = Machine(to_int, X[ftr];
            parallel=false, verbosity=verbosity - 1)
        end
        j +=1
    end

    return DiscretizerScheme(features, transformer_machines, is_ordinal)

end

function transform(transformer::Discretizer, scheme, X)

    issubset(Set(scheme.features), Set(names(X))) ||
        error("Provided DataFrame with incompatible features.")

    n_features = length(scheme.features)
    A = Array{Int}(undef, size(X, 1), n_features)

    for j in 1:n_features
        A[:,j] = transform(scheme.transformer_machines[j],
                           X[:,scheme.features[j]])
    end

    return NominalOrdinalIntArray(A, scheme.features, scheme.is_ordinal)

end

end # module


text:
"""
    VehicleState
A default type to represent an agent physical state (position, velocity).
It contains the position in the global frame, Frenet frame and the longitudinal velocity
# constructors
    VehicleState(posG::VecSE2{Float64}, v::Float64)
    VehicleState(posG::VecSE2{Float64}, roadway::Roadway, v::Float64)
    VehicleState(posG::VecSE2{Float64}, lane::Lane, roadway::Roadway, v::Float64)
    VehicleState(posF::Frenet, roadway::Roadway, v::Float64)

# fields
- `posG::VecSE2{Float64}` global position
- `posF::Frenet` lane relative position
- `v::Float64` longitudinal velocity
"""
struct VehicleState
    posG::VecSE2{Float64} # global
    posF::Frenet # lane-relative frame
    v::Float64
end

VehicleState() = VehicleState(VecSE2(), NULL_FRENET, NaN)
VehicleState(posG::VecSE2{Float64}, v::Float64) = VehicleState(posG, NULL_FRENET, v)
VehicleState(posG::VecSE2{Float64}, roadway::Roadway, v::Float64) = VehicleState(posG, Frenet(posG, roadway), v)
VehicleState(posG::VecSE2{Float64}, lane::Lane, roadway::Roadway, v::Float64) = VehicleState(posG, Frenet(posG, lane, roadway), v)
VehicleState(posF::Frenet, roadway::Roadway, v::Float64) = VehicleState(posg(posF, roadway), posF, v)

posf(veh::VehicleState) = veh.posF
posg(veh::VehicleState) = veh.posG
vel(veh::VehicleState) = veh.v
velf(veh::VehicleState) = (s=veh.v*cos(veh.posF.ϕ), t=veh.v*sin(veh.posF.ϕ))
velg(veh::VehicleState) = (x=veh.v*cos(veh.posG.θ), y=veh.v*sin(veh.posG.θ))

posf(veh::Entity) = posf(veh.state)
posg(veh::Entity) = posg(veh.state)
vel(veh::Entity)  = vel(veh.state)
velf(veh::Entity) = velf(veh.state)
velg(veh::Entity) = velg(veh.state)

Base.show(io::IO, s::VehicleState) = print(io, "VehicleState(", s.posG, ", ", s.posF, ", ", @sprintf("%.3f", s.v), ")")

function Base.write(io::IO, ::MIME"text/plain", s::VehicleState)
    @printf(io, "%.16e %.16e %.16e", s.posG.x, s.posG.y, s.posG.θ)
    @printf(io, " %d %.16e %d %d", s.posF.roadind.ind.i, s.posF.roadind.ind.t, s.posF.roadind.tag.segment, s.posF.roadind.tag.lane)
    @printf(io, " %.16e %.16e %.16e", s.posF.s, s.posF.t, s.posF.ϕ)
    @printf(io, " %.16e", s.v)
end
function Base.read(io::IO, ::MIME"text/plain", ::Type{VehicleState})
    tokens = split(strip(readline(io)), ' ')
    i = 0
    posG = VecSE2(parse(Float64, tokens[i+=1]), parse(Float64, tokens[i+=1]), parse(Float64, tokens[i+=1]))
    roadind = RoadIndex(CurveIndex(parse(Int, tokens[i+=1]), parse(Float64, tokens[i+=1])),
                        LaneTag(parse(Int, tokens[i+=1]), parse(Int, tokens[i+=1])))
    posF = Frenet(roadind, parse(Float64, tokens[i+=1]), parse(Float64, tokens[i+=1]), parse(Float64, tokens[i+=1]))
    v = parse(Float64, tokens[i+=1])
    return VehicleState(posG, posF, v)
end

"""
    Vec.lerp(a::VehicleState, b::VehicleState, t::Float64, roadway::Roadway)
Perform linear interpolation of the two vehicle states. Returns a VehicleState.
"""
function Vec.lerp(a::VehicleState, b::VehicleState, t::Float64, roadway::Roadway)
    posG = lerp(a.posG, b.posG, t)
    v = lerp(a.v, b.v, t)
    VehicleState(posG, roadway, v)
end

"""
    get_vel_s(s::VehicleState)
returns the longitudinal velocity (along the lane)
"""
get_vel_s(s::VehicleState) = s.v * cos(s.posF.ϕ) # velocity along the lane
"""
    get_vel_t(s::VehicleState)
returns the lateral velocity (⟂ to lane)
"""
get_vel_t(s::VehicleState) = s.v * sin(s.posF.ϕ) # velocity ⟂ to lane

"""
    move_along(vehstate::VehicleState, roadway::Roadway, Δs::Float64;
    ϕ₂::Float64=vehstate.posF.ϕ, t₂::Float64=vehstate.posF.t, v₂::Float64=vehstate.v)

returns a vehicle state after moving vehstate of a length Δs along its lane.
"""
function move_along(vehstate::VehicleState, roadway::Roadway, Δs::Float64;
    ϕ₂::Float64=posf(vehstate).ϕ, t₂::Float64=posf(vehstate).t, v₂::Float64=vel(vehstate)
    )

    roadind = move_along(posf(vehstate).roadind, roadway, Δs)
    try
        footpoint = roadway[roadind]
    catch
        println(roadind)
    end
    footpoint = roadway[roadind]
    posG = convert(VecE2, footpoint.pos) + polar(t₂, footpoint.pos.θ + π/2)
    posG = VecSE2(posG.x, posG.y, footpoint.pos.θ + ϕ₂)
    VehicleState(posG, roadway, v₂)
end

"""
    Vehicle
A specific instance of the Entity type defined in Records to represent Vehicles with
state `VehicleState` , definition `VehicleDef` and id `Int64`
"""
const Vehicle = Entity{VehicleState,VehicleDef,Int64}

# XXX Should this go in features
"""
    get_center(veh::Entity{VehicleState, D, I})
returns the position of the center of the vehicle
"""
get_center(veh::Entity{VehicleState, D, I}) where {D, I} = veh.state.posG
"""
    get_footpoint(veh::Entity{VehicleState, D, I})
returns the position of the footpoint of the vehicle
"""
get_footpoint(veh::Entity{VehicleState, D, I}) where {D, I} = veh.state.posG + polar(veh.state.posF.t, veh.state.posG.θ-veh.state.posF.ϕ-π/2)

"""
    get_front(veh::Entity{VehicleState, VehicleDef, I})
returns the position of the front of the vehicle
"""
get_front(veh::Entity{VehicleState, D, I}) where {D<:AbstractAgentDefinition, I} = veh.state.posG + polar(length(veh.def)/2, veh.state.posG.θ)

"""
    get_rear(veh::Entity{VehicleState, VehicleDef, I})
returns the position of the rear of the vehicle
"""
get_rear(veh::Entity{VehicleState, D, I}) where {D<:AbstractAgentDefinition, I} = veh.state.posG - polar(length(veh.def)/2, veh.state.posG.θ)

"""
    get_lane(roadway::Roadway, vehicle::Entity{S, D, I})
    get_lane(roadway::Roadway, vehicle::VehicleState)
return the lane where `vehicle` is in.
"""
function get_lane(roadway::Roadway, vehicle::Entity{S, D, I}) where {S, D, I}
    get_lane(roadway, vehicle.state)
end
function get_lane(roadway::Roadway, vehicle::VehicleState)
    lane_tag = vehicle.posF.roadind.tag
    return roadway[lane_tag]
end


text:
using JSON

function getConfig()::Dict{String, String}
    return JSON.parsefile("config.json")
end

max_stars_repo_name: tosone/julia-spider

max_issues_repo_name: tosone/julia-spider

max_forks_repo_name: tosone/julia-spider


text:
# Adapted from Distributions.jl

## Wishart

using StatsFuns: logtwo, logmvgamma

struct TuringWishart{T<:Real, ST <: Cholesky} <: ContinuousMatrixDistribution
    df::T     # degree of freedom
    chol::ST  # the Cholesky of scale matrix
    logc0::T  # the logarithm of normalizing constant in pdf
end

#### Constructors

function TuringWishart(d::Wishart)
    return TuringWishart(d.df, getchol(d.S), d.logc0)
end
getchol(p::AbstractMatrix) = cholesky(p)
getchol(p::PDMat) = p.chol
getchol(p::PDiagMat) = Cholesky(Diagonal(sqrt.(p.diag)), 'U', 0)
getchol(p::ScalMat) = Cholesky(Diagonal(fill(sqrt(p.value), p.dim)), 'U', 0)

function TuringWishart(df::Real, S::AbstractMatrix)
    p = size(S, 1)
    df > p - 1 || error("dpf should be greater than dim - 1.")
    C = getchol(S)
    return TuringWishart(df, C)
end

function TuringWishart(df::T, C::Cholesky) where {T <: Real}
    p = size(C, 1)
    df > p - 1 || error("dpf should be greater than dim - 1.")
    logc0 = _wishart_logc0(df, C)
    R = Base.promote_eltype(T, logc0)
    return TuringWishart(R(df), C, R(logc0))
end

function _wishart_logc0(df::Real, C::Cholesky)
    h_df = df / 2
    p = size(C, 1)
    -h_df * (logdet(C) + p * float(typeof(df))(logtwo)) - logmvgamma(p, h_df)
end

#### Properties

Distributions.insupport(::Type{TuringWishart}, X::Matrix) = isposdef(X)
Distributions.insupport(d::TuringWishart, X::Matrix) = size(X) == size(d) && isposdef(X)

Distributions.dim(d::TuringWishart) = size(d.chol, 1)
Base.size(d::TuringWishart) = (p = Distributions.dim(d); (p, p))
Base.size(d::TuringWishart, i) = size(d)[i]
LinearAlgebra.rank(d::TuringWishart) = Distributions.dim(d)

#### Statistics

Distributions.mean(d::TuringWishart) = d.df * Matrix(d.chol)

function Distributions.mode(d::TuringWishart)
    r = d.df - Distributions.dim(d) - 1.0
    if r > 0.0
        return Matrix(d.chol) * r
    else
        error("mode is only defined when df > p + 1")
    end
end

function Distributions.meanlogdet(d::TuringWishart)
    p = Distributions.dim(d)
    df = d.df
    v = logdet(d.chol) + p * logtwo
    for i = 1:p
        v += digamma(0.5 * (df - (i - 1)))
    end
    return v
end

function Distributions.entropy(d::TuringWishart)
    p = Distributions.dim(d)
    df = d.df
    return -d.logc0 - 0.5 * (df - p - 1) * Distributions.meanlogdet(d) + 0.5 * df * p
end

#  Gupta/Nagar (1999) Theorem 3.3.15.i
function Distributions.cov(d::TuringWishart, i::Integer, j::Integer, k::Integer, l::Integer)
    S = Matrix(d.chol)
    d.df * (S[i, k] * S[j, l] + S[i, l] * S[j, k])
end

function Distributions.var(d::TuringWishart, i::Integer, j::Integer)
    S = Matrix(d.chol)
    d.df * (S[i, i] * S[j, j] + S[i, j] ^ 2)
end

#### Evaluation

function Distributions._logpdf(d::TuringWishart, X::AbstractMatrix{<:Real})
    df = d.df
    p = Distributions.dim(d)
    Xcf = cholesky(X)
    return ((df - (p + 1)) * logdet(Xcf) - tr(d.chol \ X)) / 2 + d.logc0
end

function Distributions.logpdf(d::TuringWishart, X::AbstractArray{<:AbstractMatrix{<:Real}})
    return map(x -> logpdf(d, x), X)
end
function Distributions.logpdf(d::TuringWishart, X::AbstractArray{<:Matrix{<:Real}})
    return map(x -> logpdf(d, x), X)
end

#### Sampling
function Distributions._rand!(rng::AbstractRNG, d::TuringWishart, A::AbstractMatrix)
    Distributions._wishart_genA!(rng, A, d.df)
    unwhiten!(d.chol, A)
    A .= A * A'
end

function unwhiten!(C::Cholesky, x::StridedVecOrMat)
    cf = C.U
    lmul!(transpose(cf), x)
end

## InverseWishart

struct TuringInverseWishart{T<:Real, ST<:AbstractMatrix} <: ContinuousMatrixDistribution
    df::T     # degree of freedom
    S::ST     # Scale matrix
    logc0::T  # log of normalizing constant
end

#### Constructors

function TuringInverseWishart(d::InverseWishart)
    d = TuringInverseWishart(d.df, getmatrix(d.Ψ), d.logc0)
end
getmatrix(p::AbstractMatrix) = p
getmatrix(p::PDMat) = p.mat
getmatrix(p::PDiagMat) = Diagonal(p.diag)
getmatrix(p::ScalMat) = Diagonal(fill(p.value, p.dim))

function TuringInverseWishart(df::T, x::AbstractMatrix) where T<:Real
    Ψ = getmatrix(x)
    p = size(Ψ, 1)
    df > p - 1 || error("df should be greater than dim - 1.")
    C = getchol(x)
    logc0 = _invwishart_logc0(df, C)
    R = Base.promote_eltype(T, logc0)
    return TuringInverseWishart(R(df), Ψ, R(logc0))
end
function _invwishart_logc0(df::Real, C::Cholesky)
    h_df = df / 2
    p = size(C, 1)
    -h_df * (p * float(typeof(df))(logtwo) - logdet(C)) - logmvgamma(p, h_df)
end

#### Properties

Distributions.insupport(::Type{TuringInverseWishart}, X::Matrix) = isposdef(X)
Distributions.insupport(d::TuringInverseWishart, X::Matrix) = size(X) == size(d) && isposdef(X)

Distributions.dim(d::TuringInverseWishart) = size(d.S, 1)
Base.size(d::TuringInverseWishart) = (p = Distributions.dim(d); (p, p))
Base.size(d::TuringInverseWishart, i) = size(d)[i]
LinearAlgebra.rank(d::TuringInverseWishart) = Distributions.dim(d)

#### Statistics

function Distributions.mean(d::TuringInverseWishart)
    df = d.df
    p = Distributions.dim(d)
    r = df - (p + 1)
    if r > 0.0
        return d.S * (1.0 / r)
    else
        error("mean only defined for df > p + 1")
    end
end

Distributions.mode(d::TuringInverseWishart) = d.S * inv(d.df + Distributions.dim(d) + 1.0)

#  https://en.wikipedia.org/wiki/Inverse-Wishart_distribution#Moments
function Distributions.cov(d::TuringInverseWishart, i::Integer, j::Integer, k::Integer, l::Integer)
    p, ν, Ψ = (Distributions.dim(d), d.df, d.S)
    ν > p + 3 || throw(ArgumentError("cov only defined for df > dim + 3"))
    inv((ν - p)*(ν - p - 3)*(ν - p - 1)^2)*(2Ψ[i,j]*Ψ[k,l] + (ν-p-1)*(Ψ[i,k]*Ψ[j,l] + Ψ[i,l]*Ψ[k,j]))
end

function Distributions.var(d::TuringInverseWishart, i::Integer, j::Integer)
    p, ν, Ψ = (Distributions.dim(d), d.df, d.S)
    ν > p + 3 || throw(ArgumentError("var only defined for df > dim + 3"))
    inv((ν - p)*(ν - p - 3)*(ν - p - 1)^2)*((ν - p + 1)*Ψ[i,j]^2 + (ν - p - 1)*Ψ[i,i]*Ψ[j,j])
end

#### Evaluation

function Distributions._logpdf(d::TuringInverseWishart, X::AbstractMatrix{<:Real})
    p = Distributions.dim(d)
    df = d.df
    Xcf = cholesky(X)
    # we use the fact: tr(Ψ * inv(X)) = tr(inv(X) * Ψ) = tr(X \ Ψ)
    Ψ = d.S
    return -((df + p + 1) * logdet(Xcf) + tr(Xcf \ Ψ)) / 2 + d.logc0
end

function Distributions.logpdf(d::TuringInverseWishart, X::AbstractArray{<:AbstractMatrix{<:Real}})
    return map(x -> logpdf(d, x), X)
end
function Distributions.logpdf(d::TuringInverseWishart, X::AbstractArray{<:Matrix{<:Real}})
    return map(x -> logpdf(d, x), X)
end

#### Sampling

function Distributions._rand!(rng::AbstractRNG, d::TuringInverseWishart, A::AbstractMatrix)
    X = Distributions._rand!(rng, TuringWishart(d.df, inv(cholesky(d.S))), A)
    A .= inv(cholesky!(X))
end


text:
using Leadfield
using Winston

Xcords, Ycords, Zcords = readSRF("/Users/rluke/Documents/GITs/Leadfield.jl/data/Default50Brain.srf");

p = plotSRF(Xcords, Ycords, Zcords);
file(p, "Brain.eps", width=1200, height=600)

Xcords, Ycords, Zcords = readSRF("/Users/rluke/Documents/GITs/Leadfield.jl/data/Default50Skin.srf");

p = plotSRF(Xcords, Ycords, Zcords);
file(p, "Skin.eps", width=1200, height=600)

Xcords, Ycords, Zcords = readSRF("/Users/rluke/Documents/GITs/Leadfield.jl/data/FEM_Default50MRI.srf");

p = plotSRF(Xcords, Ycords, Zcords);
file(p, "UK.eps", width=1200, height=600)


text:
####
#### Output writer utilities
####

ext(fw::AbstractOutputWriter) = throw("Extension for $(typeof(fw)) is not implemented.")

# When saving stuff to disk like a JLD2 file, `saveproperty!` is used, which
# converts Julia objects to language-agnostic objects.
saveproperty!(file, location, p::Number)        = file[location] = p
saveproperty!(file, location, p::AbstractRange) = file[location] = collect(p)
saveproperty!(file, location, p::AbstractArray) = file[location] = Array(p)
saveproperty!(file, location, p::AbstractField) = file[location] = Array(p.data.parent)
saveproperty!(file, location, p::Function) = @warn "Cannot save Function property into $location"
saveproperty!(file, location, p) = [saveproperty!(file, location * "/$subp", getproperty(p, subp))
                                        for subp in propertynames(p)]

# Special saveproperty! so boundary conditions are easily readable outside julia.
function saveproperty!(file, location, cbcs::CoordinateBoundaryConditions)
    for endpoint in propertynames(cbcs)
        endpoint_bc = getproperty(cbcs, endpoint)
        if isa(endpoint_bc.condition, Function)
            @warn "$field.$coord.$endpoint boundary is of type Function and cannot be saved to disk!"
            file["boundary_conditions/$field/$coord/$endpoint/type"] = string(bctype(endpoint_bc))
            file["boundary_conditions/$field/$coord/$endpoint/condition"] = missing
        else
            file["boundary_conditions/$field/$coord/$endpoint/type"] = string(bctype(endpoint_bc))
            file["boundary_conditions/$field/$coord/$endpoint/condition"] = endpoint_bc.condition
        end
    end
end

# Special savepropety! for AB2 time stepper struct used by the checkpointer so
# it only saves the fields and not the tendency BCs or χ value (as they can be
# constructed by the `Model` constructor).
function saveproperty!(file, location, ts::AdamsBashforthTimestepper)
    saveproperty!(file, location * "/Gⁿ", ts.Gⁿ)
    saveproperty!(file, location * "/G⁻", ts.G⁻)
end

saveproperties!(file, structure, ps) = [saveproperty!(file, "$p", getproperty(structure, p)) for p in ps]

# When checkpointing, `serializeproperty!` is used, which serializes objects
# unless they need to be converted (basically CuArrays only).
serializeproperty!(file, location, p) = file[location] = p
serializeproperty!(file, location, p::Union{AbstractArray, AbstractField}) = saveproperty!(file, location, p)
serializeproperty!(file, location, p::Function) = @warn "Cannot serialize Function property into $location"

serializeproperties!(file, structure, ps) = [serializeproperty!(file, "$p", getproperty(structure, p)) for p in ps]

# Don't check arrays because we don't need that noise.
has_reference(T, ::AbstractArray{<:Number}) = false

# These two conditions are true, but should not necessary.
has_reference(::Type{Function}, ::AbstractField) = false
has_reference(::Type{T}, ::NTuple{N, <:T}) where {N, T} = true

"""
    has_reference(has_type, obj)

Check (or attempt to check) if `obj` contains, somewhere among its
subfields and subfields of fields, a reference to an object of type
`has_type`. This function doesn't always work.
"""
function has_reference(has_type, obj)
    if typeof(obj) <: has_type
        return true
    elseif applicable(iterate, obj) && length(obj) > 1
        return any([has_reference(has_type, elem) for elem in obj])
    elseif applicable(propertynames, obj) && length(propertynames(obj)) > 0
        return any([has_reference(has_type, getproperty(obj, p)) for p in propertynames(obj)])
    else
        return typeof(obj) <: has_type
    end
end

####
####  JLD2 output writer
####

"""
    JLD2OutputWriter{F, I, O, IF, IN, KW} <: AbstractOutputWriter

An output writer for writing to JLD2 files.
"""
mutable struct JLD2OutputWriter{F, I, O, IF, IN, KW} <: AbstractOutputWriter
        filepath :: String
         outputs :: O
        interval :: I
       frequency :: F
            init :: IF
       including :: IN
        previous :: Float64
            part :: Int
    max_filesize :: Float64
           async :: Bool
           force :: Bool
         verbose :: Bool
         jld2_kw :: KW
end

noinit(args...) = nothing

"""
    JLD2OutputWriter(model, outputs; interval=nothing, frequency=nothing, dir=".",
                     prefix="", init=noinit, including=[:grid, :coriolis, :buoyancy, :closure],
                     part=1, max_filesize=Inf, force=false, async=false, verbose=false)

Construct a `JLD2OutputWriter` that writes `label, func` pairs in `outputs` (which can be a `Dict` or `NamedTuple`)
to a JLD2 file, where `label` is a symbol that labels the output and `func` is a function of the form `func(model)`
that returns the data to be saved.

Keyword arguments
=================
- `frequency::Int`   : Save output every `n` model iterations.
- `interval::Int`    : Save output every `t` units of model clock time.
- `dir::String`      : Directory to save output to. Default: "." (current working directory).
- `prefix::String`   : Descriptive filename prefixed to all output files. Default: "".
- `init::Function`   : A function of the form `init(file, model)` that runs when a JLD2 output file is initialized.
                       Default: `noinit(args...) = nothing`.
- `including::Array` : List of model properties to save with every file. By default, the grid, equation of state,
                       Coriolis parameters, buoyancy parameters, and turbulence closure parameters are saved.
- `part::Int`        : The starting part number used if `max_filesize` is finite. Default: 1.
- `max_filesize::Int`: The writer will stop writing to the output file once the file size exceeds `max_filesize`, and
                       write to a new one with a consistent naming scheme ending in `part1`, `part2`, etc. Defaults to
                       `Inf`.
- `force::Bool`      : Remove existing files if their filenames conflict. Default: `false`.
- `async::Bool`      : Write output asynchronously. Default: `false`.
- `verbose::Bool`    : Log what the output writer is doing with statistics on compute/write times and file sizes.
                       Default: `false`.
- `jld2_kw::Dict`    : Dict of kwargs to be passed to `jldopen` when data is written.
"""
function JLD2OutputWriter(model, outputs; interval=nothing, frequency=nothing, dir=".", prefix="",
                          init=noinit, including=[:grid, :coriolis, :buoyancy, :closure],
                          part=1, max_filesize=Inf, force=false, async=false, verbose=false,
                          jld2_kw=Dict{Symbol, Any}())

    validate_interval(interval, frequency)

    mkpath(dir)
    filepath = joinpath(dir, prefix * ".jld2")
    force && isfile(filepath) && rm(filepath, force=true)

    jldopen(filepath, "a+"; jld2_kw...) do file
        init(file, model)
        saveproperties!(file, model, including)
    end

    return JLD2OutputWriter(filepath, outputs, interval, frequency, init, including,
                            0.0, part, max_filesize, async, force, verbose, jld2_kw)
end

function write_output(model, fw::JLD2OutputWriter)
    verbose = fw.verbose
    verbose && @info @sprintf("Calculating JLD2 output %s...", keys(fw.outputs))
    tc = Base.@elapsed data = Dict((name, f(model)) for (name, f)
                                    in zip(keys(fw.outputs), values(fw.outputs)))
    verbose && @info "Calculation time: $(prettytime(tc))"

    iter = model.clock.iteration
    time = model.clock.time

    filesize(fw.filepath) >= fw.max_filesize && start_next_file(model, fw)

    path = fw.filepath
    verbose && @info "Writing JLD2 output $(keys(fw.outputs)) to $path..."
    t0, sz = time_ns(), filesize(path)

    if fw.async
        @async remotecall(jld2output!, 2, path, iter, time, data, fw.jld2_kw)
    else
        jld2output!(path, iter, time, data, fw.jld2_kw)
    end

    t1, newsz = time_ns(), filesize(path)

    verbose && @info @sprintf("Writing done: time=%s, size=%s, Δsize=%s",
                              prettytime((t1-t0)/1e9), pretty_filesize(newsz), pretty_filesize(newsz-sz))

    return nothing
end

"""
    jld2output!(path, iter, time, data, kwargs)

Write the (name, value) pairs in `data`, including the simulation
`time`, to the JLD2 file at `path` in the `timeseries` group,
stamping them with `iter` and using `kwargs` when opening
the JLD2 file.
"""
function jld2output!(path, iter, time, data, kwargs)
    jldopen(path, "r+"; kwargs...) do file
        file["timeseries/t/$iter"] = time
        for (name, datum) in data
            file["timeseries/$name/$iter"] = datum
        end
    end
    return nothing
end

function start_next_file(model::Model, fw::JLD2OutputWriter)
    verbose = fw.verbose
    sz = filesize(fw.filepath)
    verbose && @info begin
        "Filesize $(pretty_filesize(sz)) has exceeded maximum file size $(pretty_filesize(fw.max_filesize))."
    end

    if fw.part == 1
        part1_path = replace(fw.filepath, r".jld2$" => "_part1.jld2")
        verbose && @info "Renaming first part: $(fw.filepath) -> $part1_path"
        mv(fw.filepath, part1_path, force=fw.force)
        fw.filepath = part1_path
    end

    fw.part += 1
    fw.filepath = replace(fw.filepath, r"part\d+.jld2$" => "part" * string(fw.part) * ".jld2")
    fw.force && isfile(fw.filepath) && rm(fw.filepath, force=true)
    verbose && @info "Now writing to: $(fw.filepath)"

    jldopen(fw.filepath, "a+"; fw.jld2_kw...) do file
        fw.init(file, model)
        saveproperties!(file, model, fw.including)
    end
end

"""
    FieldOutput([return_type=Array], field)

Returns a `FieldOutput` type intended for use with the `JLD2OutputWriter`.
Calling `FieldOutput(model)` returns `return_type(field.data.parent)`.
"""
struct FieldOutput{O, F}
    return_type :: O
          field :: F
end

FieldOutput(field) = FieldOutput(Array, field) # default
(fo::FieldOutput)(model) = fo.return_type(fo.field.data.parent)

"""
    FieldOutputs(fields)

Returns a dictionary of `FieldOutput` objects with key, value
pairs corresponding to each name and value in the `NamedTuple` `fields`.
Intended for use with `JLD2OutputWriter`.

Examples
========

```julia
julia> output_writer = JLD2OutputWriter(model, FieldOutputs(model.velocities), frequency=1);

julia> keys(output_writer.outputs)
Base.KeySet for a Dict{Symbol,FieldOutput{UnionAll,F} where F} with 3 entries. Keys:
  :w
  :v
  :u
```
"""
function FieldOutputs(fields)
    names = propertynames(fields)
    nfields = length(fields)
    return Dict((names[i], FieldOutput(fields[i])) for i in 1:nfields)
end

####
#### NetCDF output writer
####

"""
    netcdf_spatial_dimensions(::field)

Returns the dimensions associated with a field.

Examples
========

julia> netcdf_spatial_dimensions(model.velocities.u)
("xF", "yC", "zC")

julia> netcdf_spatial_dimensions(model.tracers.T)
("xC", "yC", "zC")
"""
function netcdf_spatial_dimensions(field::Field{LX, LY, LZ}) where {LX, LY, LZ}
    xdim(LX), ydim(LY), zdim(LZ)
end

xdim(::Type{Face}) = "xF"
xdim(::Type{Cell}) = "xC"
ydim(::Type{Face}) = "yF"
ydim(::Type{Cell}) = "yC"
zdim(::Type{Face}) = "zF"
zdim(::Type{Cell}) = "zC"

# This function converts an integer to a range, and nothing to a Colon
get_slice(n::Integer) = n:n
get_slice(n::UnitRange) = n
get_slice(n::Nothing) = Colon()

"""
    write_grid(model; filename="grid.nc", mode="c",
               compression=0, attributes=Dict(), slice_kw...)

Writes a grid.nc file that contains all the dimensions of the domain.

Keyword arguments
=================

    - `filename::String`  : File name to be saved under
    - `mode::String`      : Netcdf file is opened in either clobber ("c") or append ("a") mode (Default: "c" )
    - `compression::Int`  : Defines the compression level of data (0-9, default 0)
    - `attributes::Dict`  : Attributes (default: Dict())
"""
function write_grid(model; filename="./grid.nc", mode="c",
                    compression=0, attributes=Dict(), slice_kw...)

    dimensions = Dict(
        "xC" => collect(model.grid.xC),
        "yC" => collect(model.grid.yC),
        "zC" => collect(model.grid.zC),
        "xF" => collect(model.grid.xF)[1:end-1],
        "yF" => collect(model.grid.yF)[1:end-1],
        "zF" => collect(model.grid.zF)[1:end-1]
    )

    dim_attrib = Dict(
        "xC" => ["longname" => "Locations of the cell centers in the x-direction.", "units" => "m"],
        "yC" => ["longname" => "Locations of the cell centers in the y-direction.", "units" => "m"],
        "zC" => ["longname" => "Locations of the cell centers in the z-direction.", "units" => "m"],
        "xF" => ["longname" => "Locations of the cell faces in the x-direction.",   "units" => "m"],
        "yF" => ["longname" => "Locations of the cell faces in the y-direction.",   "units" => "m"],
        "zF" => ["longname" => "Locations of the cell faces in the z-direction.",   "units" => "m"]
    )

    # Applies slices to the dimensions d
    for (d, slice) in slice_kw
        if String(d) in keys(dimensions)
            dimensions[String(d)] = dimensions[String(d)][get_slice(slice)]
        end
    end

    # Writes the sliced dimensions to the specified netcdf file
    Dataset(filename, mode, attrib=attributes) do ds
        for (dimname, dimarray) in dimensions
            defVar(ds, dimname, dimarray, (dimname,),
                   compression=compression, attrib=dim_attrib[dimname])
        end
    end

    return nothing
end

"""
    NetCDFOutputWriter <: AbstractOutputWriter

An output writer for writing to NetCDF files.
"""
mutable struct NetCDFOutputWriter <: AbstractOutputWriter
             filename :: String
              dataset :: Any
              outputs :: Dict
             interval :: Union{Nothing, AbstractFloat}
            frequency :: Union{Nothing, Int}
              clobber :: Bool
               slices :: Dict
   len_time_dimension :: Int
             previous :: Float64
end

"""
    NetCDFOutputWriter(model, outputs; interval=nothing, frequency=nothing, filename=".",
                                   clobber=true, global_attributes=Dict(), output_attributes=nothing, slice_kw...)

Construct a `NetCDFOutputWriter` that writes `label, field` pairs in `outputs` (which can be a
`Dict` or `NamedTuple`) to a NC file, where `label` is a symbol that labels the output and
`field` is a field from the model (e.g. `model.velocities.u`).

Keyword arguments
=================

    - `filename::String`         : Directory to save output to. Default: "." (current working directory).
    - `frequency::Int`           : Save output every `n` model iterations.
    - `interval::Int`            : Save output every `t` units of model clock time.
    - `clobber::Bool`            : Remove existing files if their filenames conflict. Default: `false`.
    - `compression::Int`         : Determines the compression level of data (0-9, default 0)
    - `global_attributes::Dict`  : Dict of model properties to save with every file (deafult: Dict())
    - `output_attributes::Dict`  : Dict of attributes to be saved with each field variable (reasonable
                                   defaults are provided for velocities, temperature, and salinity)
    - `slice_kw`                 : `dimname = Union{OrdinalRange, Integer}` will slice the dimension `dimname`.
                                   All other keywords are ignored.
                                   e.g. `xC = 3:10` will only produce output along the dimension `xC` between
                                   indices 3 and 10 for all fields with `xC` as one of their dimensions.
                                   `xC = 1` is treated like `xC = 1:1`.
                                   Multiple dimensions can be sliced in one call. Not providing slices writes
                                   output over the entire domain.
"""

function NetCDFOutputWriter(model, outputs; interval=nothing, frequency=nothing, filename=".",
                        clobber=true, global_attributes=Dict(), output_attributes=nothing, compression=0, slice_kw...)

    validate_interval(interval, frequency)

    mode = clobber ? "c" : "a"

    # Initiates the output file with dimensions
    write_grid(model; filename=filename, mode=mode,
               compression=compression, attrib=global_attributes, slice_kw...)

    # Opens the same output file for writing fields from the user-supplied variable outputs
    dataset = Dataset(filename, "a")

    # Creates an unliimited dimension "Time"
    defDim(dataset, "Time", Inf); sync(dataset)
    defVar(dataset, "Time", Float64, ("Time",)); sync(dataset)
    len_time_dimension = 0 # Number of time-steps so far

    if isa(output_attributes, Nothing)
        output_attributes = Dict("u" => Dict("longname" => "Velocity in the x-direction", "units" => "m/s"),
                            "v" => Dict("longname" => "Velocity in the y-direction", "units" => "m/s"),
                            "w" => Dict("longname" => "Velocity in the z-direction", "units" => "m/s"),
                            "T" => Dict("longname" => "Temperature", "units" => "K"),
                            "S" => Dict("longname" => "Salinity", "units" => "g/kg"))
    end

    # Initiates empty Float32 arrays for fields from the user-supplied variable outputs
    for (fieldname, field) in outputs
        dtype = eltype(interiorparent(field))
        defVar(dataset, fieldname, dtype, (netcdf_spatial_dimensions(field)...,"Time"),
               compression=compression, attrib=output_attributes[fieldname])
    end
    sync(dataset)

    # Stores slices for the dimensions of each output field
    slices = Dict{String, Vector{Union{OrdinalRange,Colon}}}()
    for (fieldname, field) in outputs
        slices[fieldname] = slice(field; slice_kw...)
    end

    return NetCDFOutputWriter(filename, dataset, outputs, interval,
                          frequency, clobber, slices, len_time_dimension, 0.0)
end

# Closes the outputwriter
Base.close(fw::NetCDFOutputWriter) = close(fw.dataset)

function read_output(fw::NetCDFOutputWriter, fieldname)
    ds = Dataset(fw.filename,"r")
    field = ds[fieldname][:,:,:,end]
    close(ds)
    return field
end

"""
    slice(field; slice_kw...)

For internal use only. Returns a slice for a field based on its dimensions and the supplied slices in `slice_kw`.
"""
function slice(field; slice_kw...)
    slice_array = Vector{Union{AbstractRange,Colon}}()
    for dim in netcdf_spatial_dimensions(field)
        slice = get_slice(haskey(slice_kw, Symbol(dim)) ? slice_kw[Symbol(dim)] : nothing)
        push!(slice_array, slice)
    end
    return slice_array
end

"""
    write_output(model, OutputWriter)

For internal user only. Writes output to the netcdf file at specified intervals. Increments the `Time` dimension every time an output is written to the file.
"""
function write_output(model, fw::NetCDFOutputWriter)
    fw.len_time_dimension += 1
    fw.dataset["Time"][fw.len_time_dimension] = model.clock.time
    for (fieldname, field) in fw.outputs
        fw.dataset[fieldname][:, :, :, fw.len_time_dimension] = view(interiorparent(field), fw.slices[fieldname]...)
    end
    sync(fw.dataset)
    return nothing
end

####
#### Checkpointer
####

"""
    Checkpointer{I, T, P, A} <: AbstractOutputWriter

An output writer for checkpointing models to a JLD2 file from which models can be restored.
"""
mutable struct Checkpointer{I, T, P, A} <: AbstractOutputWriter
         frequency :: I
          interval :: T
          previous :: Float64
               dir :: String
            prefix :: String
        properties :: P
    has_array_refs :: A
             force :: Bool
           verbose :: Bool
end

"""
    Checkpointer(model; frequency=nothing, interval=nothing, dir=".", prefix="checkpoint",
                        force=false, verbose=false,
                        properties = [:architecture, :boundary_conditions, :grid, :clock,
                                      :eos, :constants, :closure, :velocities, :tracers,
                                      :timestepper])

Construct a `Checkpointer` that checkpoints the model to a JLD2 file every so often as
specified by `frequency` or `interval`. The `model.clock.iteration` is included in the
filename to distinguish between multiple checkpoint files.

Note that extra model `properties` can be safely specified, but removing crucial properties
such as `:velocities` will make restoring from the checkpoint impossible.

The checkpoint file is generated by serializing model properties to JLD2. However,
functions cannot be serialized to disk (at least not with JLD2). So if a model property
contains a reference somewhere in its hierarchy it will not be included in the checkpoint
file (and you will have to manually restore them).

Keyword arguments
=================
- `frequency::Int`   : Save output every `n` model iterations.
- `interval::Int`    : Save output every `t` units of model clock time.
- `dir::String`      : Directory to save output to. Default: "." (current working directory).
- `prefix::String`   : Descriptive filename prefixed to all output files. Default: "checkpoint".
- `force::Bool`      : Remove existing files if their filenames conflict. Default: `false`.
- `verbose::Bool`    : Log what the output writer is doing with statistics on compute/write times and file sizes.
                       Default: `false`.
- `properties::Array`: List of model properties to checkpoint.
"""
function Checkpointer(model; frequency=nothing, interval=nothing, dir=".", prefix="checkpoint", force=false,
                      verbose=false, properties = [:architecture, :boundary_conditions, :grid, :clock, :coriolis,
                                                   :buoyancy, :closure, :velocities, :tracers, :timestepper])

    validate_interval(frequency, interval)

    # Grid needs to be checkpointed for restoring to work.
    :grid ∉ properties && push!(properties, :grid)

    has_array_refs = Symbol[]

    for p in properties
        isa(p, Symbol) || @error "Property $p to be checkpointed must be a Symbol."
        p ∉ propertynames(model) && @error "Cannot checkpoint $p, it is not a model property!"

        if has_reference(Function, getproperty(model, p))
            @warn "model.$p contains a function somewhere in its hierarchy and will not be checkpointed."
            filter!(e -> e != p, properties)
        end

        has_reference(AbstractField, getproperty(model, p)) && push!(has_array_refs, p)
    end

    mkpath(dir)
    return Checkpointer(frequency, interval, 0.0, dir, prefix, properties, has_array_refs, force, verbose)
end

function write_output(model, c::Checkpointer)
    filepath = joinpath(c.dir, c.prefix * string(model.clock.iteration) * ".jld2")
    c.verbose && @info "Checkpointing to file $filepath..."

    t0 = time_ns()
    jldopen(filepath, "w") do file
        file["checkpointed_properties"] = c.properties
        file["has_array_refs"] = c.has_array_refs

        serializeproperties!(file, model, filter(e -> !(e ∈ c.has_array_refs), c.properties))
        saveproperties!(file, model, filter(e -> e ∈ c.has_array_refs, c.properties))
    end

    t1, sz = time_ns(), filesize(filepath)
    c.verbose && @info "Checkpointing done: time=$(prettytime((t1-t0)/1e9)), size=$(pretty_filesize(sz))"
end

defaultname(::Checkpointer, nelems) = :checkpointer
_arr(::CPU, a) = a
_arr(::GPU, a) = CuArray(a)

function restore_fields!(model, file, arch, fieldset; location="$fieldset")
    if fieldset == :timestepper
        restore_fields!(model.timestepper, file, arch, :Gⁿ; location="timestepper/Gⁿ")
        restore_fields!(model.timestepper, file, arch, :G⁻; location="timestepper/G⁻")
    else
        for p in propertynames(getproperty(model, fieldset))
            getproperty(getproperty(model, fieldset), p).data.parent .= _arr(arch, file[location * "/$p"])
        end
    end
end

"""
    restore_from_checkpoint(filepath; kwargs=Dict())

Restore a model from the checkpoint file stored at `filepath`. `kwargs` can be passed
to the `Model` constructor, which can be especially useful if you need to manually
restore forcing functions or boundary conditions that rely on functions.
"""
function restore_from_checkpoint(filepath; kwargs=Dict())
    file = jldopen(filepath, "r")

    cps = file["checkpointed_properties"]
    has_array_refs = file["has_array_refs"]

    # Restore model properties that were just serialized.
    # We skip fields that contain structs and restore them after model creation.
    for p in cps
        if p ∉ has_array_refs
            kwargs[p] = file["$p"]
        end
    end

    model = Model(; kwargs...)

    # Now restore fields.
    for p in cps
        if p in has_array_refs
            restore_fields!(model, file, model.architecture, p)
        end
    end

    close(file)

    return model
end


text:
#ДАНО: Робот - в произвольной клетке ограниченного прямоугольного поля с перегородками
#РЕЗУЛЬТАТ: Робот - в исходном положении, и все клетки поля промакированы
include("library.jl")
function zakras_with_obstacles(r::Robot)
    sides = moveAndReturnDirections!(r)
    marks!(r)
    moveToStartplace!(r)
    moveToBeginplace!(r, sides)
end

max_stars_repo_name: listova-nastya/mirea-progs

max_issues_repo_name: listova-nastya/mirea-progs

max_forks_repo_name: listova-nastya/mirea-progs


text:
using PooksoftAlphaVantageDataStore

# ok, so lets build by user object -
my_current_dir = pwd()   # where am I?
path_to_config_file = my_current_dir*"/test/configuration/Configuration.json"

# build the api user model -
user_model = build_api_user_model(path_to_config_file) |> checkresult

# make an API call -
ticker_symbol = "ally"
data_type = :json
df = execute_sts_quote_api_call(user_model, ticker_symbol; data_type = data_type, logger=nothing) |> checkresult


text:
# small offset along ray direction to prevent ray from reintersecting the surface it just intersected. Optical pathlength is measured in mm so this is much smaller than a wavelength of light.
const RAY_OFFSET = 1e-6
# if the power of a ray drops below this value return nothing and do not continue the ray tracing recursion
const POWER_THRESHOLD = 1e-6
# if the ray is traced recursively more than this many times then the trace will be terminated and the ray ignored (mostly happens with TIR)
const TRACE_RECURSION_LIMIT = 256

# triangulated surfaces are imprecise mostly due to floating point precision, particularly if using trig functions
# expanding every triangle by a fraction resolves missed intersections
const TRIANGULATION_EXPANSION = 2 * eps(Float64)
# we also need to extend the bounds of the surface to ensure we don't miss anythgin at the edges of the surface
const TRIANGULATION_EXTENSION = 0.02

# for accelerated surface the triangulated surface can often be very imprecise, this will mean rays which start close
# to the surface might miss intersection because the ray origin lies above the triangulated surface but below the true surface
# (where for this example the ray direction is straight upwards), offsetting the ray fixes this, and as we check the alpha
# of the intersection point afterwards there are no problems
const ACCEL_SURF_RAY_OFFSET = 0.1

# settings for visualization
const MIN_VIS_TRI_SIZE = 1e-3
const MESH_PRECISION = 1e-3
const VIS_RECURSION_LIMIT = 64

# number of M/N to precompute for various variables involved in the QType calculations
const QTYPE_PRECOMP = 16

# maximum number of orders for a thin grating interface
const GRATING_MAX_ORDERS = 10

# LensAssemblyN methods will be generated at compile time (faster on first run) for N up to this number
const PREGENERATED_LENS_ASSEMBLY_SIZE = 10


text:
export Monoflop

function Monoflop(name::Symbol, duration::H; start_delay=zero(H), M=Int, add_reset=false) where H
    P = Symbol[]
    T = Symbol[]
    m₀= M[]
    eft=H[]
    lft=H[]
    trans_probs=Float64[]

    R =  Tuple{Int,Int,M}[]
    ΔF = Tuple{Int,Int,M}[]

    function add_place!(name,initial)
        push!(P, Symbol(name))
        push!(m₀, M(initial))
        length(P)
    end

    function add_transition!(name,e,l=e,p=1.0)
        push!(T, Symbol(name))
        push!(eft, H(e))
        push!(lft, H(l))
        push!(trans_probs, p)
        length(T)
    end

    # add places
    p_on,p_off = add_place!.((("$(name)_$(state)") for state in (on,:off)),(0,1))

    # add transitions and arcs
    t_start,t_stop = add_transition!.((("$(name)_$(trans)") for trans in (:start,:stop)),(start_delay,duration))
    push!(ΔF,
        (p_off, t_start, -1),
        (p_on,  t_start,  1),
        (p_off, t_stop,   1),
        (p_on,  t_stop,  -1)
    )

    # add reset-related transitions and arcs
    if add_reset
        # add transitions
        t_reset = add_transition!("$(name)_reset",zero(H),zero(H),1.0)
        # add state-changing arcs
        push!(ΔF,
            (p_on, t_reset, -1),
            (p_off, t_reset, 1)
        )
    end

    ΔF = sparse(collect.(zip(ΔF...))..., length(P), length(T))
    R  = isempty(R) ? zero(ΔF) : sparse(collect.(zip(R...))..., length(P), length(T))

    TPN(P,T,R,ΔF,eft,lft,trans_probs,m₀)
end


text:
module LocalREPLHistory

HISTORY_FILE = ".julia_history"
GITIGNORE = ".gitignore"

function init()
    open(joinpath(homedir(), ".julia/config/startup.jl"), "a") do fp
        print(fp, """
using LocalREPLHistory; LocalREPLHistory.use_localhistoryfile()
""")
    end
end

function use_localhistoryfile(targetdir = pwd())
    if HISTORY_FILE in readdir(targetdir)
        ENV["JULIA_HISTORY"] = joinpath(pwd(), HISTORY_FILE)
    end
end

function add_localhistoryfile(targetdir = pwd(); append_gitignore = true)
    checkfiles = filter(x -> isfile(x), readdir(targetdir))
    if !(HISTORY_FILE in checkfiles)
        history_file = joinpath(targetdir, HISTORY_FILE)
        open(history_file, "w") do fp
            write(fp, "")
        end
    end
    if append_gitignore && GITIGNORE in checkfiles
        gifile = joinpath(targetdir, GITIGNORE)
        gilines = readlines(gifile,keep=true)
        if !(HISTORY_FILE in gilines)
            open(gifile, "a") do fp
                if length(gilines) != 0 && gilines[end][end] != '\n'
                    print(fp, "\n")
                end
                print(fp, HISTORY_FILE * "\n")
            end
        end
    end
    use_localhistoryfile(targetdir)
end

add = add_localhistoryfile

function rm_localhistoryfile(targetdir = pwd())
    checkfiles = filter(x -> isfile(x), readdir(targetdir))
    if HISTORY_FILE in checkfiles
        Base.Filesystem.rm(joinpath(targetdir, HISTORY_FILE))
    end
    if GITIGNORE in checkfiles
        gifile = joinpath(targetdir, GITIGNORE)
        none_historyfilegi = filter(x -> x != HISTORY_FILE, readlines(gifile))
        if length(none_historyfilegi) != 0
            open(gifile, "w") do fp
                for l in none_historyfilegi
                    println(fp, l)
                end
            end
        end
    end
end

rm = rm_localhistoryfile

end # module


text:
function makeplot(func,picdir)
    PyPlot.clf()
    f=getfield(Main,Symbol(func))
    gridplot(f(), Plotter=PyPlot)
    fname=joinpath(picdir,func*".svg")
    PyPlot.savefig(fname)
    isfile(fname)
end

function makeplots(picdir)

    makeplot("interval_from_vector",picdir)
    makeplot("interval_localref",picdir)
    makeplot("interval_multiregion",picdir)
    makeplot("interval_subgrid",picdir)
    makeplot("rectangle",picdir)
    makeplot("rectangle_localref",picdir)
    makeplot("rectangle_multiregion",picdir)
    makeplot("rectangle_subgrid",picdir)
    makeplot("quadrilateral",picdir)

end