text:
### A Pluto.jl notebook ###
# v0.14.5

using Markdown
using InteractiveUtils

# This Pluto notebook uses @bind for interactivity. When running this notebook outside of Pluto, the following 'mock version' of @bind gives bound variables a default value (instead of an error).
macro bind(def, element)
    quote
        local el = $(esc(element))
        global $(esc(def)) = Core.applicable(Base.get, el) ? Base.get(el) : missing
        el
    end
end

# ╔═╡ 97e18b69-4eb6-42a9-beef-04697156a7ab
using PlutoUI, JuMP

# ╔═╡ 2f845732-c095-11eb-14e8-57175fb5038e
lines = readlines(open("./benchmark_problems/Class B/INST13.txt"))

# ╔═╡ 12f3987f-ea86-46a2-852f-65ce5886cbe2
split_line(line) = filter(x->x!="", split(line, " "))

# ╔═╡ 2c521bfc-683b-4337-b62e-a5586d3d06cb
function parse_int(i)
	try
		parse(Int, i)
	catch ArgumentError
		parse(Int, replace(i, ".00"=>""))
	end
end

# ╔═╡ d4175f5d-d4e9-4862-9c60-83a74af1dc1c
parse_line(line) = map(n->parse_int(n), split_line(line))

# ╔═╡ 0707ad63-038b-4b0f-9c48-310502701915
file = map(parse_line, lines[2:end])

# ╔═╡ dca8d05b-f878-4fb6-bb8e-cef62dd2249a
m_blocks, n_rows, n_groups = file[1]

# ╔═╡ f1d972de-aa4a-4e5e-91d6-890bd4b21fcd
groups_bounds_row = file[2]

# ╔═╡ 75894d1c-b1c4-4050-8b52-19d837e2bc31
data = file[3:end]

# ╔═╡ cf5f3e19-7087-4e67-bb3e-d02e2b33efeb
block_start_index(i) = (i-1)*(n_rows+1) + 1

# ╔═╡ dc9221b9-09b8-48e9-b46b-8cad6cb0b778
@bind row Slider(1:n_groups)

# ╔═╡ 39ba33bb-7113-4f35-9f19-676178e210cb
n_groups

# ╔═╡ c9fe0f02-8f6a-4584-88ca-d0640fab9be2
n_i = n_rows

# ╔═╡ 2ef7e159-9bff-4b36-92a3-52c43db1c8fe
m = n_groups

# ╔═╡ d7205749-58ab-4831-8797-ac3d7a175f37
falses(m)

# ╔═╡ c06e268b-09a2-482f-98fa-6ecbd8ee9fe1
n = m_blocks # size of set of disjoint groups

# ╔═╡ 9c4c419c-0517-4ae9-9a48-7368444085fe
function load_block(i, n_groups, data)
	d = data[block_start_index(i):block_start_index(i+1)-1]
	@assert first(first(d)) == i
	items = d[2:end]
	#exactly one of each element in this group must be turned on
	#in order to represent this in MPS form, we need to one-hot encode the class index
	index = falses(n_groups)
	index[i] = true

	#then combine the one-hot encoding with the resource requirements
	resource_reqs = map(x->vcat(x[2:end], index), items)

	#strip the profit vector, then return the (profit, resources) tuple
	profit = map(x->first(x), items)
	profit, resource_reqs
end

# ╔═╡ 1c10debb-ee04-45e2-81fc-9fe4d9362bc7
current_block = load_block(row, n_groups, data)

# ╔═╡ e764cd4e-fc74-4561-8265-e303b2892c88
length(current_block)

# ╔═╡ 8652e037-b90b-4044-a618-730f1f554b2e
function make_model(profit, resource_reqs, available_res)
	m = Model()
	x = @variable(m, x[1:length(profit)])
	add_resource_constraint(m, x, r, a) =
		@constraint(m, sum(x .* r) <= a)
	map((i, row)->add_resource_constraint(m, x, row, a[i]), resource_reqs)
	@objective(m, Max, sum(x .* profit))
	m
end

# ╔═╡ d3938ef5-e05a-47c3-b4c6-bc966245acd8
blocks = map(i->load_block(i, n_groups, data), 1:n_groups)

# ╔═╡ c35b78db-1c40-4313-a580-83f48680a7e1
profits = vcat(map(b->first(b), blocks)...)

# ╔═╡ 1db24586-4b11-4b75-bfee-1015644527de
resource_reqs = vcat(map(b->b[2], blocks)...)

# ╔═╡ 8fd748aa-6aff-4dcc-a5fd-c50d52295dc5
resource_limits = vcat(groups_bounds_row..., ones(n_groups))

# ╔═╡ Cell order:
# ╠═2f845732-c095-11eb-14e8-57175fb5038e
# ╠═12f3987f-ea86-46a2-852f-65ce5886cbe2
# ╠═d4175f5d-d4e9-4862-9c60-83a74af1dc1c
# ╠═2c521bfc-683b-4337-b62e-a5586d3d06cb
# ╠═0707ad63-038b-4b0f-9c48-310502701915
# ╠═dca8d05b-f878-4fb6-bb8e-cef62dd2249a
# ╠═f1d972de-aa4a-4e5e-91d6-890bd4b21fcd
# ╠═75894d1c-b1c4-4050-8b52-19d837e2bc31
# ╠═cf5f3e19-7087-4e67-bb3e-d02e2b33efeb
# ╠═1c10debb-ee04-45e2-81fc-9fe4d9362bc7
# ╠═d7205749-58ab-4831-8797-ac3d7a175f37
# ╠═dc9221b9-09b8-48e9-b46b-8cad6cb0b778
# ╠═39ba33bb-7113-4f35-9f19-676178e210cb
# ╠═e764cd4e-fc74-4561-8265-e303b2892c88
# ╠═c9fe0f02-8f6a-4584-88ca-d0640fab9be2
# ╠═2ef7e159-9bff-4b36-92a3-52c43db1c8fe
# ╠═c06e268b-09a2-482f-98fa-6ecbd8ee9fe1
# ╠═9c4c419c-0517-4ae9-9a48-7368444085fe
# ╠═8652e037-b90b-4044-a618-730f1f554b2e
# ╠═d3938ef5-e05a-47c3-b4c6-bc966245acd8
# ╠═c35b78db-1c40-4313-a580-83f48680a7e1
# ╠═1db24586-4b11-4b75-bfee-1015644527de
# ╠═8fd748aa-6aff-4dcc-a5fd-c50d52295dc5
# ╠═97e18b69-4eb6-42a9-beef-04697156a7ab

max_stars_repo_name: Dibillilia/JuliaOptHeuristics

max_issues_repo_name: Dibillilia/JuliaOptHeuristics

max_forks_repo_name: Dibillilia/JuliaOptHeuristics


text:
import FileIO
import ImageMagick
import Images

function _bytes_to_image(file_extension::AbstractString, bytes::AbstractVector{UInt8})
    original_directory = pwd()
    temp_directory = mktempdir()
    atexit(() -> rm(temp_directory; force = true, recursive = true))
    temp_image_file = joinpath(temp_directory, "image.$(file_extension)")
    rm(temp_image_file; force = true, recursive = true)
    open(temp_image_file, "w") do io
        write(io, bytes)
    end
    image = FileIO.load(temp_image_file)
    cd(original_directory)
    rm(temp_image_file; force = true, recursive = true)
    rm(temp_directory; force = true, recursive = true)
    return image
end

ext: jl


text:
### A Pluto.jl notebook ###
# v0.19.8

using Markdown
using InteractiveUtils

# This Pluto notebook uses @bind for interactivity. When running this notebook outside of Pluto, the following 'mock version' of @bind gives bound variables a default value (instead of an error).
macro bind(def, element)
    quote
        local iv = try Base.loaded_modules[Base.PkgId(Base.UUID("6e696c72-6542-2067-7265-42206c756150"), "AbstractPlutoDingetjes")].Bonds.initial_value catch; b -> missing; end
        local el = $(esc(element))
        global $(esc(def)) = Core.applicable(Base.get, el) ? Base.get(el) : iv(el)
        el
    end
end

# ╔═╡ 4712d6fd-efb3-4f46-9160-c82797a3b5e4
# ╠═╡ show_logs = false
begin
	using Logging
	global_logger(NullLogger())
	display("")
end

# ╔═╡ f0545c67-5cfd-438f-a9ef-92c35ebaefa4
# ╠═╡ show_logs = false
#Set-up packages
begin

	using DataFrames, HTTP, CSV, Dates, Plots, PlutoUI, Printf, LaTeXStrings, HypertextLiteral, ShortCodes

	gr();
	Plots.GRBackend()

	#Define html elements
	nbsp = html"&nbsp" #non-breaking space
	vspace = html"""

"""
	br = html"
"

	#Sets the width of cells, caps the cell width by 90% of screen width
	#(setting overwritten by cell below)
	# @bind screenWidth @htl("""
	#

	#
	#

	# """)

	# cellWidth= min(1000, screenWidth*0.9)
	# @htl("""
	#
	# """)

	#Sets the width of the cells
	#begin
	#	html"""


text:
using DiffEqCallbacks, OrdinaryDiffEq, Test

# Non-negative ODE examples
#
# Reference:
# Shampine, L.F., S. Thompson, J.A. Kierzenka, and G.D. Byrne,
# "Non-negative solutions of ODEs," Applied Mathematics and Computation Vol. 170, 2005,
# pp. 556-569.
# https://www.mathworks.com/help/matlab/math/nonnegative-ode-solution.html

"""
Absolute value function

```math
\\frac{du}{dt} = -|u|
```
with initial condition ``u₀=1``, and solution

```math
u(t) = u₀*e^{-t}
```
for positive initial values ``u₀``.
"""
function absval(du,u,p,t)
    du[1] = -abs(u[1])
end
analytic(u₀, p, t) = u₀*exp(-t)
ff = ODEFunction(absval,analytic=analytic)
prob_absval = ODEProblem(ff, [1.0], (0.0, 40.0))

# naive approach leads to large errors
naive_sol_absval = solve(prob_absval, BS3())
@test naive_sol_absval.errors[:l∞] > 9e4
@test naive_sol_absval.errors[:l2] > 1.3e4

# general domain approach
# can only guarantee approximately non-negative values
function g(resid,u)
    resid[1] = u[1] < 0 ? -u[1] : 0
end
general_sol_absval = solve(prob_absval, BS3(); callback=GeneralDomain(g, [1.0]))
@test_broken all(x -> x[1] ≥ -10*eps(), general_sol_absval.u)
@test general_sol_absval.errors[:l∞] < 9.9e-5
@test general_sol_absval.errors[:l2] < 4.3e-5
@test general_sol_absval.errors[:final] < 4.3e-18

# test non-autonomous function
g_t(resid, u, p, t) = g(resid, u)

general_t_sol_absval = solve(prob_absval, BS3(); callback=GeneralDomain(g_t, [1.0]))
@test general_sol_absval.t ≈ general_t_sol_absval.t &&
    general_sol_absval.u ≈ general_t_sol_absval.u

# positive domain approach
# can guarantee non-negative values
positive_sol_absval = solve(prob_absval, BS3(); callback=PositiveDomain([1.0]))
@test_broken all(x -> x[1] ≥ 0, positive_sol_absval.u)
@test_broken general_sol_absval.t ≈ positive_sol_absval.t
@test general_sol_absval.errors[:l∞] ≈ positive_sol_absval.errors[:l∞]

# specify abstol as array or scalar
positive_sol_absval2 = solve(prob_absval, BS3(); callback=PositiveDomain([1.0], abstol=[1e-6]))
@test positive_sol_absval.t ≈ positive_sol_absval2.t &&
    positive_sol_absval.u ≈ positive_sol_absval2.u
positive_sol_absval3 = solve(prob_absval, BS3(); callback=PositiveDomain([1.0], abstol=1e-6))
@test positive_sol_absval.t ≈ positive_sol_absval3.t && positive_sol_absval.u ≈ positive_sol_absval3.u

# specify scalefactor
positive_sol_absval4 = solve(prob_absval, BS3(); callback=PositiveDomain([1.0], scalefactor=0.2))
@test length(positive_sol_absval.t) < length(positive_sol_absval4.t)
@test positive_sol_absval.errors[:l2] > positive_sol_absval4.errors[:l2]

"""
Knee problem

```math
\\frac{du}{dt} = \epsilon^{-1}(1-t-u)u
```

with initial condition ``u0=1``, and generally ``0 < \epsilon << 1``.
Here ``\epsilon=1e-6``. Then the solution approaches ``u=1-t`` for ``t<1``
and ``u=0`` for ``t>1``.
"""
function knee(du,u,p,t)
    du[1] = 1e6*(1-t-u[1])*u[1]
end

prob_knee = ODEProblem(knee, [1.0], (0.0, 2.0))

# unfortunately callbacks do not work with solver CVODE_BDF which is comparable to ode15s
# used in MATLAB example, so we use Rodas5
naive_sol_knee = solve(prob_knee, Rodas5())
@test naive_sol_knee[1, end] ≈ -1.0 atol=1e-5

# positive domain approach
positive_sol_knee = solve(prob_knee, Rodas5(); callback=PositiveDomain([1.0]))
@test_broken all(x -> x[1] ≥ 0, positive_sol_knee.u)
@test positive_sol_knee[1, end] ≈ 0.0

## Now test on out-of-place equations
r, K = 1.1, 10.0
logistic(u,p,t) = u*r*(1-u/K)
t = (0.0, 20.0)
logistic_p = ODEProblem(logistic, 0.02, t)
logistic_s = solve(logistic_p, Tsit5())
logistic_s_positive = solve(logistic_p, Tsit5(), callback=PositiveDomain())


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

name = "UCX"
version = v"1.7.0"

# Collection of sources required to complete build
sources = [
    ArchiveSource("https://github.com/openucx/ucx/releases/download/v1.7.0/ucx-1.7.0.tar.gz",
                  "6ab81ee187bfd554fe7e549da93a11bfac420df87d99ee61ffab7bb19bdd3371"),
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir/ucx-*
./configure --prefix=${prefix} \
    --build=${MACHTYPE} \
    --host=${target} \
    --disable-logging \
    --disable-debug \
    --disable-assertions \
    --disable-params-check \
    --disable-static
# For a bug in `src/uct/sm/cma/Makefile` that I did't have the time to look
# into, we have to build with `V=1`
make -j${nproc} V=1
make install
"""

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

# The products that we will ensure are always built
products = [
    LibraryProduct("libuct", :libuct),
    LibraryProduct("libucm", :libucm),
    LibraryProduct("libucp", :libucp),
    LibraryProduct("libucs", :libucs),
    ExecutableProduct("ucx_info", :ucx_info),
    ExecutableProduct("ucx_perftest", :ucx_perftest),
    ExecutableProduct("ucx_read_profile", :ucx_read_profile),
]

# Dependencies that must be installed before this package can be built
dependencies = [
    Dependency(PackageSpec(name="CompilerSupportLibraries_jll", uuid="e66e0078-7015-5450-92f7-15fbd957f2ae")),
    Dependency(PackageSpec(name="NUMA_jll", uuid="7f51dc2b-bb24-59f8-b771-bb1490e4195d")),
]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies)


text:
function shercliff(;nx::Int=3, ny::Int=3, Re::Float64 = 10.0, Ha::Float64 = 10.0,
    f_u = nothing, resultsfile = nothing)

  L = 1.0
  N = Ha^2/Re
  K = Ha / (1-0.825*Ha^(-1/2)-Ha^(-1))
  ∂p∂z = -Re * K / L^3

  if f_u === nothing
    _f_u = VectorValue(0.0,0.0, -∂p∂z) * L/Re^2 # Assumed ν=L=1 by default
  else
    _f_u = f_u
  end
  g_u = VectorValue(0.0,0.0,0.0)
  g_j = VectorValue(0.0,0.0,0.0)
  B = VectorValue(0.0,1.0,0.0)

  # Discretization
  order = 2
  domain = (-1.0,1.0,-1.0,1.0,0.0,0.1)
  map(x) = VectorValue(sign(x[1])*(abs(x[1])*0.5)^0.5,
                       sign(x[2])*(abs(x[2])*0.5)^0.5, x[3]*sqrt(2)/2)*2/sqrt(2)

  dirichlet_tags_u = append!(collect(1:20),[23,24,25,26])
  dirichlet_tags_j = append!(collect(1:20),[23,24,25,26])

  partition=(nx,ny,3)
  model = CartesianDiscreteModel(domain,partition;isperiodic=(false,false,true))

  labels = get_face_labeling(model)
  add_tag_from_tags!(labels,"dirichlet_u",dirichlet_tags_u)
  add_tag_from_tags!(labels,"dirichlet_j",dirichlet_tags_j)

  Vu = FESpace(model, ReferenceFE(lagrangian,VectorValue{3,Float64},order);
      conformity=:H1, dirichlet_tags="dirichlet_u")

  Vp = FESpace(model, ReferenceFE(lagrangian,Float64,order-1,space=:P);
      conformity=:L2, constraint=:zeromean)

  Vj = FESpace(model, ReferenceFE(raviart_thomas,Float64,order-1);
      conformity=:Hdiv, dirichlet_tags="dirichlet_j")

  Vφ = FESpace(model, ReferenceFE(lagrangian,Float64,order-1,space=:Q);
      conformity=:L2, constraint=:zeromean)

  U = TrialFESpace(Vu,g_u)
  P = TrialFESpace(Vp)
  J = TrialFESpace(Vj,g_j)
  Φ = TrialFESpace(Vφ)

  Y = MultiFieldFESpace([Vu, Vp, Vj, Vφ])
  X = MultiFieldFESpace([U, P, J, Φ])

  trian = Triangulation(model)
  degree = 2*(order)
  dΩ = Measure(trian,degree)

  res(x,y) = InductionlessMHD.dimensionless_residual(x, y, Re, N, B, _f_u,dΩ)
  jac(x,dx,y) = InductionlessMHD.dimensionless_jacobian(x, dx, y, Re, N, B,dΩ)

  op  = FEOperator(res,jac,X,Y)

  nls = NLSolver(;
    show_trace=true, method=:newton, linesearch=BackTracking())
  solver = FESolver(nls)

  xh = solve(solver,op)

  if resultsfile != nothing
    uh, ph, jh, φh = xh
    writevtk(trian, resultsfile,
      cellfields=["uh"=>uh, "ph"=>ph, "jh"=>jh, "φh"=>φh])
  end

  (xh, trian, dΩ)
end

function analytical_shercliff_u(a::Float64,       # semi-length of side walls
                                b::Float64,       # semi-length of Hartmann walls
                                t_w::Float64,     # wall thickness
                                σ_w::Float64,     # wall conductivity
                                σ::Float64,       # fluid conductivity
                                μ::Float64,       # fluid viscosity
                                grad_pz::Float64, # presure gradient
                                Ha::Float64,      # Hartmann number
                                n::Int,           # number of sumands included in Fourier series
                                x)                # evaluation point
  l = b/a
  ξ = x[1]/a
  η = x[2]/a

  d_B = t_w*σ_w/(a*σ)

  V = 0.0; V0=0.0;
  for k in 0:n
    α_k = (k + 0.5)*π/l
    N = (Ha^2 + 4*α_k^2)^(0.5)
    r1_k = 0.5*( Ha + N)
    r2_k = 0.5*(-Ha + N)

    num1 = d_B*r2_k + (1-exp(-2*r2_k))/(1+exp(-2*r2_k))
    num2 = (exp(-r1_k*(1-η))+exp(-r1_k*(1+η)))/2.0
    den1 = d_B*N *(1+exp(-2*r1_k))/2.0
    den2 = (1-exp(-2*(r1_k+r2_k)))/(1+exp(-2*r2_k))
    V2 = (num1 * num2)/(den1 + den2)

    num1 = d_B*r1_k + (1-exp(-2*r1_k))/(1+exp(-2*r1_k))
    num2 = (exp(-r2_k*(1-η)) + exp(-r2_k*(1+η)))/2.0
    den1 = d_B*N*(1+exp(-2*r2_k))/2
    den2 = (1-exp(-2(r1_k+r2_k)))/(1+exp(-2*r1_k))
    V3 = (num1 * num2)/(den1 + den2)

    V += 2*(-1)^k*cos(α_k * ξ)/(l*α_k^3) * (1-V2-V3)
  end
  u_z = V/μ * (-grad_pz) * a^2

  return VectorValue(0.0*u_z,0.0*u_z,u_z)
end

function analytical_shercliff_j(a::Float64,       # semi-length of side walls
                                b::Float64,       # semi-length of Hartmann walls
                                t_w::Float64,     # wall thickness
                                σ_w::Float64,     # wall conductivity
                                σ::Float64,       # fluid conductivity
                                μ::Float64,       # fluid viscosity
                                grad_pz::Float64, # presure gradient
                                Ha::Float64,      # Hartmann number
                                n::Int,           # number of sumands included in Fourier series
                                x)                # evaluation point
  l = b/a
  ξ = x[1]/a
  η = x[2]/a

  d_B = t_w*σ_w/(a*σ)

  H_dx = 0.0; H_dy = 0.0
  for k in 0:n
    α_k = (k + 0.5)*π/l
    N = sqrt(Ha^2 + 4*α_k^2)
    r1_k = 0.5*( Ha + N)
    r2_k = 0.5*(-Ha + N)

    num1 = d_B*r2_k + (1-exp(-2*r2_k))/(1+exp(-2*r2_k))
    num2 = (exp(-r1_k*(1-η)) - exp(-r1_k*(1+η)))/2.0
    num2_dy = (exp(-r1_k*(1-η))*(r1_k/a) + exp(-r1_k*(1+η))*(r1_k/a))/2.0
    den1 = d_B*N*(1+exp(-2*r1_k))/2
    den2 = (1-exp(-2*(r1_k+r2_k)))/(1+exp(-2*r2_k))
    H2 = (num1 * num2)/(den1 + den2)
    H2_dy = (num1 * num2_dy)/(den1 + den2)

    num1 = d_B*r1_k + (1-exp(-2*r1_k))/(1+exp(-2*r1_k))
    num2 = (exp(-r2_k*(1-η)) - exp(-r2_k*(1+η)))/2.0
    num2_dy = (exp(-r2_k*(1-η))*(r2_k/a) + exp(-r2_k*(1+η))*(r2_k/a))/2.0
    den1 = d_B*N*(1+exp(-2*r2_k))/2
    den2 = (1-exp(-2*(r1_k+r2_k)))/(1+exp(-2*r1_k))
    H3 = (num1 * num2)/(den1 + den2)
    H3_dy = (num1 * num2_dy)/(den1 + den2)

    H_dx += -2*(-1)^k * sin(α_k * ξ)/(a*l*α_k^2) * (H2 - H3)
    H_dy += 2*(-1)^k * cos(α_k * ξ)/(l*α_k^3) * (H2_dy - H3_dy)
  end
  j_x = a^2*σ^0.5 / μ^0.5 * (-grad_pz) * H_dy
  j_y = a^2*σ^0.5 / μ^0.5 * (-grad_pz) * (-H_dx)

  return VectorValue(j_x,j_y,0.0)
end


text:
# ------------------------------------------------------------------
# Licensed under the MIT License. See LICENSE in the project root.
# ------------------------------------------------------------------

# helper function for table view
function viewtable(table, rows, cols)
  t = Tables.columns(table)
  v = map(cols) do c
    col = Tables.getcolumn(t, c)
    c => view(col, rows)
  end
  (; v...)
end

"""
    DataView(sdata, inds, vars)

Return a view of spatial data `sdata` at `inds` and `vars`.
"""
struct DataView{𝒮,I,V} <: AbstractData
  data::𝒮
  inds::I
  vars::V
end

domain(dv::DataView) = view(domain(dv.data), dv.inds)
values(dv::DataView) = viewtable(values(dv.data), dv.inds, dv.vars)

Base.collect(dv::DataView) = georef(values(dv), coordinates(dv))

# -----------------------------------
# specialize methods for performance
# -----------------------------------
coordinates!(buff::AbstractVector, dv::DataView, ind::Int) =
  coordinates!(buff, dv.data, dv.inds[ind])

# ---------
# VIEW API
# ---------

Base.view(dv::DataView, inds::AbstractVector{Int}) =
  DataView(dv.data, dv.inds[inds], dv.vars)
Base.view(dv::DataView, vars::AbstractVector{Symbol}) =
  DataView(dv.data, dv.inds, vars)
Base.view(dv::DataView, inds, vars) =
  DataView(dv.data, dv.inds[inds], vars)

# ------------
# IO methods
# ------------
function Base.show(io::IO, dv::DataView)
  N = ncoords(dv)
  T = coordtype(dv)
  n = nelms(dv)
  print(io, "$n DataView{$T,$N}")
end

function Base.show(io::IO, ::MIME"text/plain", dv::DataView)
  𝒟 = domain(dv)
  𝒯 = values(dv)
  s = Tables.schema(𝒯)
  vars = zip(s.names, s.types)
  println(io, 𝒟)
  println(io, "  variables")
  varlines = ["    └─$var ($V)" for (var,V) in vars]
  print(io, join(sort(varlines), "\n"))
end


text:
function constructportfolios(;
  refreshcomp::Bool = true,
  refreshtretcrsp::Bool = true,
  refreshtvolmcapcrsp::Bool = true,
  refreshdataseries::Bool = true,
  refreshportfolios::Bool = true,
  portname::String = PORT_NAME,
  panelname::String = PANEL_NAME,
  datapath::String=DATA_PATH,
  injlsstream::Function = IN_JLS_STREAM,
  outjlsstream::Function = OUT_JLS_STREAM,
  outpath::String = OUT_PATH,
  trialrunonly::Bool = false)::DataFrame

  local panel::DataFrame
  local port::DataFrame

  #load the leverage data set
  local dailyunlesstrial::Symbol = trialrunonly ? :monthly : :daily
  local cqports::CQuartileConstructions = leverageconstructions()
  panelname = "$panelname-$dailyunlesstrial"
  portname = "$portname-$dailyunlesstrial"

  if refreshportfolios
    local panel = formdataseries(
      refreshcomp = refreshcomp,
      refreshtvolmcapcrsp = refreshtvolmcapcrsp,
      refreshtretcrsp = refreshtretcrsp,
      refreshdataseries=refreshdataseries,
      trialrunonly = trialrunonly)

    panel = makepanelmonthly(panel, crsptype=dailyunlesstrial)
    port = formfactors(panel, crsptype=dailyunlesstrial)
    cqports = merge(cqports, filteredleverageconstructions())
    #println(describe(port))
    port = mergeoutside(port)
    #println(describe(port))

    #merge and verify integrity
    println("Size of panel: $(size(panel))")
    panel = join(panel, port, on=:date, kind=:inner)
    (s->finiteormissing!(panel, s)).([:ret; :lret; T1_CONTROL10])
    #winsorizekeydata!(panel, fields2winsorize = [:ret; :lret; :ret12m; :LLret12m; :vol252d; :LLvol252d; T1_CONTROL10])

    checkdfkey(panel, [:date, :lpermno])
    println("Size of panel: $(size(panel))")

    checksum::Float64 = computechecksum(panel) + computechecksum(port)
    @info "Panel constructed with checksum $(checksum)"

    condensewithlast!(panel, condense=false)
    if TEST_OUTPUT
      d = Date(2015,12,31)
      CSV.write("output\\$(panelname)_$d.csv", panel[(cd->cd==d).(panel.date), :])
      CSV.write("output\\$(panelname)_ts_$d.csv", panel[(end - 20_000):end, :])
      CSV.write("output\\$portname.csv", port[:,
        filter(s::Symbol->string(s)[1:2] ≠ "S_", names(port))])
    end

    #WARNING These rows save space, but are not essential
    select!(port, Not(Fns(cqports)))
    retainedpanelcols = unique([:gvkey; :permno; :date; :fyear; :linkeffdate; :linkenddate;
      :adate; :begindate; :enddate; names(port); Fgroups(cqports); Fvals(cqports);
      :ret; :lret; T1_CONTROL10; :ret12m; :LLret12m; :vol252d; :Lvol252d;
      :LLvol252d; :annflag; (s->Symbol(:LD, s)).(CQUART_FVALS); CQUART_FVALS_BRS_D;])
    select!(panel, retainedpanelcols)

    sort!(panel, [:permno, :date])
    outjlsstream("$datapath\\$panelname.jls.lz4") do s
      serialize(s, panel)
    end
    #CSV.write("$datapath\\$portname.csv")
  #error("A-OK!")
  else
    injlsstream("$datapath\\$panelname.jls.lz4") do s
      panel = deserialize(s) #NOTE: DEBUG ONLY
    end
  end

  return panel
end

#=@inline function finiteormissing!(df::AbstractDataFrame, s::Symbol)::Nothing
  df[!,s] .= ((f::Union{Missing, Real})->((!ismissing(f)) && isfinite(f)) ? f : missing).(df[!,s])

  return nothing
end=#

#use this to create the default portoflio constructions
function leverageconstructions()
  local cqports::CQuartileConstructions

  local cqports = CQuartileConstructions(CQUART_FVALS, CQUART_FCONTROLS,
    Fweightons=CQUART_FWEIGHTONS, namevec=CQUART_FNAMES)
  local cqportsalt = CQuartileConstructions(CQUART_FVALS, CQUART_FCONTROLS_ALT,
    Fweightons=CQUART_FWEIGHTONS, namevec=CQUART_FNAMES_ALT)
  cqports = merge(cqports, cqportsalt)

  local cqportsd = CQuartileConstructions(CQUART_FVALS_D, CQUART_FCONTROLS_D,
    Fweightons=CQUART_FWEIGHTONS_D, namevec=CQUART_FNAMES_D)
  cqports = merge(cqports, cqportsd)
  local cqportsaltd = CQuartileConstructions(CQUART_FVALS_D, CQUART_FCONTROLS_ALT_D,
    Fweightons=CQUART_FWEIGHTONS_D, namevec=CQUART_FNAMES_ALT_D)
  cqports = merge(cqports, cqportsaltd)

  local cqportsrmd = CQuartileConstructions(CQUART_FVALS_RMD, CQUART_FCONTROLS_RMD,
    Fweightons=CQUART_FWEIGHTONS_RMD, namevec=CQUART_FNAMES_RMD)
  cqports = merge(cqports, cqportsrmd)
  local cqportsaltrmd = CQuartileConstructions(CQUART_FVALS_RMD, CQUART_FCONTROLS_ALT_RMD,
    Fweightons=CQUART_FWEIGHTONS_RMD, namevec=CQUART_FNAMES_ALT_RMD)
  cqports = merge(cqports, cqportsaltrmd)

  local cqportsbrs = CQuartileConstructions(CQUART_FVALS_BRS_D, CQUART_FCONTROLS_BRS_D,
    Fweightons=CQUART_FWEIGHTONS_BRS_D, namevec=CQUART_FNAMES_BRS_D)
  cqports = merge(cqports, cqportsbrs)

  return cqports
end

function filteredleverageconstructions()
  local cqfiltered::CQuartileConstructions

  cqfiltered = CQuartileConstructions(CQUART_FVALS_D, CQUART_FCONTROLS_D,
    Fweightons=CQUART_FWEIGHTONS_D, namevec=CQUART_FNAMES_FD)

  return cqfiltered
end

#computes a checksum of a dataframe. Used to check consistency.
function computechecksum(df::DataFrame; checksum::Float64 = 0.0)::Float64
  for c ∈ eachcol(df)
    if eltype(c) <: Union{Missing, Real}
      checksum += sum((x->ismissing(x) ? 0.0 : isfinite(x) ? x : 0.0).(c))
    end
  end
  checksum = log(checksum)

  return checksum
end

#this set of functions computes the returns from the panel
#accomplished via dispath
const NandRet = Tuple{Int, MFloat64}

function NandRet(panel::AbstractDataFrame, Fweighton::Symbol; Fret::Symbol=FRET)::NandRet
  spanel::SubDataFrame = view(
    panel, completecases(panel[!, [Fret, Fweighton]]), [Fret, Fweighton])

  N::Int = size(spanel,1)
  (N==0) && return (0,missing)
  return NandRet(N, spanel[!,Fret], spanel[!,Fweighton])
end

function NandRet(panel::AbstractDataFrame, ::Nothing; Fret::Symbol=FRET)::NandRet

  rets::T where T<:SubArray = view(panel[!, Fret], (!ismissing).(panel[!, Fret]))
  N::Int = length(rets)
  (N==0) && (return (0,missing))
  return NandRet(N, rets)
end

NandRet(N::Int, rets::AbstractVector)::NandRet = (N, sum(rets)/N)
NandRet(N::Int, rets::AbstractVector, weights::AbstractVector
  )::NandRet = (N, sum(rets .* weights)/sum(weights))

#gets a flat return file
#assumes the existance of a return file
function formfactors(panel::DataFrame; Fret::Symbol=FRET, crsptype::Symbol)::DataFrame

  local port::DataFrame
  local cqports::CQuartileConstructions = merge(leverageconstructions(), filteredleverageconstructions())

  local fgrps::Vector{Symbol} = Fgroups(cqports)
  local debugsuffix="$crsptype"

  #get the monthly returns

  #this will hold all of the portfolios
  #get the returns
  port = DataFrame(date = unique(panel.date))
  sort!(port, :date)

  Nport::Int = size(port, 1)

  #form the necessary columns
  for pc ∈ cqports
    port[!, pc.name] = missings(MFloat64, Nport)
    port[!, pc.Fn] = missings(MInt, Nport)
  end
  #index the return dates
  portrows::Dict = Dict(r.date=>r for r ∈ eachrow(port))

  #good point to check our answers
  if TEST_OUTPUT
    d = Date(2015,12,31)
    CSV.write("output\\panel$(debugsuffix)_$d-mid.csv", panel[(cd->cd==d).(panel.date), :])
  end

  #form the portfolio returns
  #print("PROFILING: portfolio returns loop - ")
  spanels::GroupedDataFrame = groupby(panel, :date)
  @mpar for i ∈ 1:length(spanels)
    spanel::SubDataFrame = spanels[i]
    r::DataFrameRow = portrows[spanel.date[1]]

    #iterate through our portfolio constructions
    for pc ∈ cqports
      r[pc.name] = 0.0
      r[pc.Fn] = 0
      for grp::Int ∈ 1:pc.Ngroups #add the contribution of each group to the factor return
        sspanel = view(spanel, (g::MInt->(!ismissing(g)) && (g == grp)).(spanel[!,pc.Fgroup]),:)

        local N::Int
        local val::MFloat64
        (N, val) = NandRet(sspanel, pc.Fweighton)
        r[pc.Fn] += N
        r[pc.name] += val * pc.groupweights[grp]
      end
    end
  end

  return port
end

function constructportfolios!(lev::DataFrame)
  #first construct the portfolios
  local cqports::CQuartileConstructions = leverageconstructions()

  N::Int = size(lev,1)

  for pc ∈ cqports
    lev[!, pc.Fgroup] = Vector{MInt}(undef, N)
    lev[!, Symbol(:c, pc.Fgroup)] = Vector{MInt}(undef, N) #this is a temporary field holding the 4 within
  end

  #assign the gorups and weights
  lastcontrol::NSymbol = nothing
  one2four::Vector{Int} = collect(1:4)

  #use this to avoid allocation within the innerloop
  sortperms::Vector{Vector{Int}} = [Vector{Int}(undef, 4) for i ∈ 1:Threads.nthreads()]
  for pc ∈ cqports

    csym::Symbol = Symbol(:c, pc.Fgroup)
    #first do the control sort
    if lastcontrol ≠ pc.Fcontrol
      sort!(lev, pc.Fcontrol)
      lastcontrol = pc.Fcontrol
    end

    #complete cases only to make the logic easier
    local essentialfields::Vector{Symbol} = filter!(
      s->!isnothing(s), [pc.Fval, pc.Fcontrol, pc.Fweighton])
    local slev::SubDataFrame = view(lev,completecases(lev[:,essentialfields]), :)
    sslevs = groupby(slev, :fyear)
    for i ∈ 1:length(sslevs) #for each fiscal year
      local sslev::SubDataFrame = sslevs[i]
      local N::Int = size(sslev,1)

      sslev[1:N, csym] .= (i->i÷4).(0:(N-1)) #create the four-group blocks
      s3levs = groupby(sslev, csym)
      @mpar for j ∈ 1:length(s3levs) #for each block
        local s3lev::SubDataFrame = s3levs[j]
        (size(s3lev,1) ≠ 4) && continue

        tid::Int = Threads.threadid()
        sortperm!(sortperms[tid], s3lev[!, pc.Fval])
        #println("vals: $(s3lev[!, pc.Fval])")
        #println("order: $(sortperms[tid])\n")
        s3lev[sortperms[tid], pc.Fgroup] .= one2four
      end
      #N>10 && error("look at the results")
    end
  end

  if TEST_OUTPUT
    y::Int = 2015
    CSV.write("output\\levtemptest_$y.csv", lev[(cd->cd==y).(lev.fyear), :])
  end
  #println("*****THIS IS IT*****:\n", describe(lev))
end

#replaces certain group entries with missing
function filterportfoliogroup!(lev::DataFrame, Fgrp::Symbol, filterfunction::Function,
  Ffiltered::Symbol = Symbol(:F, Fgrp))

  lev[!,Ffiltered] = Vector{MInt}(undef, size(lev,1))

  for (i,r) ∈ enumerate(eachrow(lev))
    r[Ffiltered] = filterfunction(r, Fgrp)
  end

  return nothing
end

const DVOL_FILTER_FIELD = :LDvol252d
#only allows for grp 1 if leverage decreases and group 4 if leverage increases
@inline function filterbyvol(r::DataFrameRow, Fgrp::Symbol, FDvol::Symbol = DVOL_FILTER_FIELD)::MInt
  grp::MInt = r[Fgrp]

  Dvol::MFloat64 = r[FDvol]
  if ismissing(Dvol) || ismissing(grp)
    return missing
  end

  #only care about groups 1 and 4
  if (grp == 2 ||
    grp == 3 ||
    (grp == 1 && Dvol < 0) ||
    (grp == 4 && Dvol > 0))

    return grp
  end

  return missing
end

function doubleincreasegroups!(lev::DataFrame;
  grpfields::Vector{Symbol} = CQUART_FGROUPS_D,
  filteredgrpfields::Vector{Symbol} = CQUART_FGROUPS_FD,
  filterfunction::Function = filterbyvol)::Nothing

  ((Fgrp::Symbol, Ffiltered::Symbol)->
    filterportfoliogroup!(lev, Fgrp, filterfunction, Ffiltered)).(grpfields, filteredgrpfields)

  return nothing
end

max_stars_repo_name: clintonTE/CCA

max_issues_repo_name: clintonTE/CCA

max_forks_repo_name: clintonTE/CCA


text:
# TODO: Merge _lasso and lasso
function _lasso(data::ModelSelection.ModelSelectionData)
    return _lasso!(data)
end

# TODO: Merge _lasso! and lasso!
function _lasso!(data::ModelSelection.ModelSelectionData)
    res = lasso!(data)
    res[1].extras[:lasso] = Dict()
    res[1].extras[:lasso][:betas] = res[2]

    return res[1]
end

function lasso(data::ModelSelection.ModelSelectionData)
    lasso!(ModelSelection.copy_data(data))
end

function lasso!(data::ModelSelection.ModelSelectionData; addextrasflag=true)
    betas, lambda = lassoselection(data)

    if isnothing(betas)
        return data, map(b -> true, data.expvars)
    end

    data.extras[:lasso_betas] = betas

    vars = map(b -> b != 0, betas)
    lassonumvars = size(filter(b -> b != 0, betas), 1)

    if data.intercept
        vars[ModelSelection.get_column_index(:_cons, data.expvars)] = true
    end

    data.expvars = data.expvars[vars]
    data.expvars_data = data.expvars_data[:,vars]

    if( addextrasflag )
        data = addextras(data, lassonumvars, betas, lambda)
    end

    return data, vars
end

function computablevars(nvars::Int)
    return 20
    min(Int(floor(log(2,Sys.total_memory()/2 ^30) + 21)), nvars)
end

function lassoselection(data)
    data = ModelSelection.filter_data_by_empty_values(data)
    data = ModelSelection.convert_data(data)
    nvars = computablevars(size(data.expvars,1))

    if nvars >= size(data.expvars,1)
        return nothing, nothing
    end

    path = glmnet(data.expvars_data, data.depvar_data; nlambda=1000)

    best = 1
    for (i, cant) in enumerate(nactive(path.betas))
        if cant >= nvars
            if cant == nvars
                best = i
            end
            break
        end
        best = i
    end

    return path.betas[:, best], path.lambda[best]
end


text:
using SolveDSGE, Plots;
if 1==1
    open("Det_INV.txt", "w") do io
        println(io,"# inv\n")
        println(io,"states:\nk \nend\n")
        println(io,"jumps:\ni \nend\n")
        println(io,"shocks:\nend\n")
        println(io,"parameters:\nδ=0.10, ρ=0.15, z=0.27 \nend\n")
        println(io,"equations:")
        println(io,"k(+1) = (1-δ)*k + i")
        println(io,"i(+1) = ((ρ+δ-z)/(1-δ)) +((1+ρ)/(1-δ))*i")
        println(io,"end")
    end;
    ###########################################################################
    path = joinpath(@__DIR__,"Det_INV.txt")
    process_model(path)
    ngm = retrieve_processed_model(path)
    ###########################################################################
    # Steady State
    ###########################################################################
    tol = 1e-8
    maxiters = 1000
    x = .25*ones(ngm.number_variables) #x = [0.25,0.25] #x = [0.25,0.25, .25] #x = [0.25,0.25, .25, .25]
    ss = compute_steady_state(ngm, x, tol, maxiters)
    ss = ss.zero
    ss ≈ [K_SS, I_SS]
    sss = ss[1:ngm.number_states] #ss for states only
    ###########################################################################
    # Perturbation, Order = 1,2,3
    # soln_fo/soln_so/soln_to
    ###########################################################################
    if 1==1
        g0 = ss;          # point about which to perturb the model (generally the steady state)
        cutoff = 1.0      # param that separates unstable from stable eigenvalues
        # order = "first" # order of the perturbation
        N   = PerturbationScheme(g0,cutoff,"first")   # Klein (2000)
        NN  = PerturbationScheme(g0,cutoff,"second")  # Gomme and Klein (2011)
        NNN = PerturbationScheme(g0,cutoff,"third")   # Binning (2013) w refinement from Levintal (2017)

        soln_fo = solve_model(ngm, N)
        soln_so = solve_model(ngm, NN)
        soln_to = solve_model(ngm, NNN)
    end

    # ProjectionSchemes: ChebyshevSchemes/SmolyakSchemes/PiecewiseLinearSchemes
    g0 = ss;                  #initial guess
    dom = [sss*2.0; sss*0.01] #domain for state vars. 2*n_s array.  #row 1 upper vals. row 2 lower vals
    tf = 1e-8;  #inner loop tol
    tv = 1e-6;  #outer loop tol
    n_thr = 4; #num threads to use.
    #nq = 9; #number quadrature points. Int. #Don't need if deterministic!!!!
    ###########################################################################
    # Projection, ChebyshevSchemeDet
    # soln_cha
    ###########################################################################
    if 1==1
        ng = chebyshev_nodes; #node generator: chebyshev_nodes/chebyshev_extrema
        #nn1 = [21]; #number of nodes for each state var #vec integers nn2 = [71];

        nn1 = 21 * ones(ngm.number_states) .|> Int64
        nn2 = 71 * ones(ngm.number_states) .|> Int64
        no1 = 6; #order of Cheby poly

        P    = ChebyshevSchemeDet(g0,ng,nn1,no1,dom,tf,tv,maxiters)
        PP   = ChebyshevSchemeDet(g0,ng,nn1,4,  dom,tf,tv,maxiters)
        PPP  = ChebyshevSchemeDet(g0,ng,nn2,6,  dom,tf,tv,maxiters)

        #
        soln_cha = solve_model(ngm,P)
        soln_cha = solve_model(ngm,P,n_thr)
        soln_chb = solve_model(ngm,PP,n_thr)
        soln_chc = solve_model(ngm,PPP,n_thr)

    end

    ###########################################################################
    # Projection, SmolyakSchemeDet            # soln_sa
    # Projection, PiecewiseLinearSchemeStoch  # soln_pl
    # Compose                                 # soln_c1
    ###########################################################################
    # Smolyak
    if 1==1
        ng = chebyshev_gauss_lobatto; #node generator: chebyshev_gauss_lobatto/clenshaw_curtis_equidistant
        nl1 = 3; #number of layers to be used in the approximation
        no1 = 6; #order of Cheby poly
        S = SmolyakSchemeDet(g0,ng,nl1,dom,tf,tv,maxiters)
        soln_sa = solve_model(ngm,S,n_thr)
    end

    # PL
    if 1==1
        nn1 = [21]; #number of nodes for each state var #vec integers
        PL = PiecewiseLinearSchemeDet(g0,nn1,dom,tf,tv,maxiters)
        soln_pl = solve_model(ngm,PL,n_thr)
    end
    # Compose
    soln_c1 = solve_model(ngm,soln_to,PPP)
    soln_c2 = solve_model(ngm,soln_c1,PP)
    # Simulation
    if 1 == 1
        T = 100 # num Sim pds
        # T=10
        ω = 1.0 - 0.2; # distance from steady-state
        #ω = 1.0 - 0.75;

        sim_d1   = simulate(soln_fo, sss*ω,  T)
        sim_d2   = simulate(soln_so, sss*ω,  T)
        sim_d3   = simulate(soln_to, sss*ω,  T)
        sim_d4   = simulate(soln_cha, sss*ω, T)
        sim_d5   = simulate(soln_sa, sss*ω,  T)
        sim_d6   = simulate(soln_pl, sss*ω,  T)
        sim_d7   = simulate(soln_c1, sss*ω,  T)

        # Closed Form Sol

        lv = reshape(ngm.variables,1,ngm.number_variables)
        p  = plot(legend=:bottomright);
        #p0 = plot!(1:T,cf,      title="Closed Form",    lab="");
        p0 = plot!([K_SS I_SS],  seriestype = :hline, lab="", color="grey", l=:dash)
        p1 = plot(p0, 1:T,sim_d1', title="Pert. 1st Order", lab=lv)
        p2 = plot(p0, 1:T,sim_d2', title="Pert. 2nd Order", lab=lv)
        p3 = plot(p0, 1:T,sim_d3', title="Pert. 3rd Order", lab=lv)
        p4 = plot(p0, 1:T,sim_d4', title="Proj. Cheby 1",   lab=lv)
        p5 = plot(p0, 1:T,sim_d5', title="Proj. Smol  1",   lab=lv)
        p6 = plot(p0, 1:T,sim_d6', title="Proj. PL    1",   lab=lv)
        p7 = plot(p0, 1:T,sim_d7', title="Proj. & Pert",    lab=lv)
        #p8 = plot!(1:T,sim_d8', title="Proj. & Pert",   lab=lv)
        plot(p1,p2,p3,p4,p5,p6, size=1.25 .* (600, 400))
    end
    #
end


text:
function load_functions(p::AbstractModel)
    if     p.numerics.jacobian == :symbolic
        jac_type = jacobian_symbolic
        load_func = load_functions_symbolic
    elseif p.numerics.jacobian == :AD
        jac_type = jacobian_AD
        load_func = load_functions_forward_diff
    end

    ## Begin loading functions
    initial_guess!, f_alg!, f_diff!, J_y!, J_y_alg!, θ_keys = load_func(p)

    append!(p.cache.θ_keys, θ_keys)
    append!(p.cache.θ_tot,  zeros(length(θ_keys)))

    funcs = model_funcs(initial_guess!, f_diff!, f_alg!, J_y!, J_y_alg!)

    return funcs
end

function get_saved_model_version(p::AbstractModel)
    """
    Gets the version of PETLION used to generate the saved sol
    """
    str = strings_directory_func(p)
    str *= "/info.txt"

    out = readline(str)
    out = replace(out, "PETLION version: v" => "")
    # convert this to a tuple
    out = (Meta.parse.(split(out,"."))...,)
end

function remove_model_files(p::AbstractModel)
    """
    Removes symbolic files for the sol
    """
    str = strings_directory_func(p) * "/"
    for x in readdir(str)
        rm(str*x)
    end
end

function load_functions_symbolic(p::AbstractModel)
    dir = strings_directory_func(p) * "/"
    files_exist = isdir(dir)

    if files_exist
        file_version = get_saved_model_version(p)

        # have there been any breaking changes since creating the functions?
        no_breaking_changes = (file_version[1] == PETLION_VERSION[1]) && (file_version[2] == PETLION_VERSION[2])

        if !no_breaking_changes
            @warn "Breaking updates encountered: re-evaluating sol..."
            remove_model_files(p)
        end

    else
        no_breaking_changes = false
    end

    if files_exist && no_breaking_changes && options[:SAVE_SYMBOLIC_FUNCTIONS]
        ## residuals
        initial_guess! = include(dir * "initial_guess.jl")
        f_alg!         = include(dir * "f_alg.jl")
        f_diff!        = include(dir * "f_diff.jl")
        J_y!_func      = include(dir * "J_y.jl")
        J_y_alg!_func  = include(dir * "J_y_alg.jl")

        ## Jacobian
        @load dir * "J_sp.jl" J_y_sp θ_keys
    else
        Y0Func,res_algFunc,res_diffFunc,jacFunc,jac_algFunc,J_y_sp,θ_keys = generate_functions_symbolic(p)

        initial_guess! = eval(Y0Func)
        f_alg!         = eval(res_algFunc)
        f_diff!        = eval(res_diffFunc)
        J_y!_func      = eval(jacFunc)
        J_y_alg!_func  = eval(jac_algFunc)
    end

    J_y!_sp     = sparse(J_y_sp...)
    J_y_alg!_sp = @inbounds J_y!_sp[p.N.diff+1:end,p.N.diff+1:end]

    J_y!     = jacobian_symbolic(J_y!_func, J_y!_sp)
    J_y_alg! = jacobian_symbolic(J_y_alg!_func, J_y_alg!_sp)

    return initial_guess!, f_alg!, f_diff!, J_y!, J_y_alg!, θ_keys
end

function generate_functions_symbolic(p::AbstractModel; verbose=options[:SAVE_SYMBOLIC_FUNCTIONS])

    if verbose println("Creating the functions for $p\n\nMay take a few minutes...") end

    θ_sym, Y_sym, YP_sym, t_sym, SOC_sym, X_applied, γ_sym, p_sym, θ_keys = get_symbolic_vars(p)

    ## Y0 function
    if verbose println("1/4: Making initial guess function") end
    Y0_sym = _symbolic_initial_guess(p_sym, SOC_sym, θ_sym, X_applied)

    ## batteryModel function
    if verbose println("2/4: Making symbolic sol") end
    res = _symbolic_residuals(p_sym, t_sym, Y_sym, YP_sym)

    θ_sym_slim, θ_keys_slim = get_only_θ_used_in_model(θ_sym, θ_keys, res, Y0_sym)

    Y0Func = build_function(Y0_sym, SOC_sym, θ_sym_slim, X_applied, fillzeros=false, checkbounds=false)[2]

    ## jacobian
    if verbose println("3/4: Making symbolic Jacobian") end
    Jac, jacFunc, J_sp = _symbolic_jacobian(p, res, t_sym, Y_sym, YP_sym, γ_sym, θ_sym, θ_sym_slim, verbose=verbose)

    if verbose println("4/4: Making initial condition functions") end
    res_algFunc, res_diffFunc = _symbolic_initial_conditions_res(p, res, t_sym, Y_sym, YP_sym, θ_sym, θ_sym_slim)

    jac_algFunc = _symbolic_initial_conditions_jac(p, Jac, t_sym, Y_sym, YP_sym, γ_sym, θ_sym_slim)

    J_y_sp = (findnz(J_sp)..., p.N.tot-1, p.N.tot)

    θ_keys = θ_keys_slim

    if options[:SAVE_SYMBOLIC_FUNCTIONS]
        dir = strings_directory_func(p; create_dir=true) * "/"

        ## file info
        str = model_info(p)

        open(dir * "info.txt", "w") do f
            write(f, str)
        end

        save_string(x) = (string(x) |> remove_comments |> rearrange_if_statements |> join)

        write(dir * "initial_guess.jl", save_string(Y0Func))
        write(dir * "J_y.jl",           save_string(jacFunc))
        write(dir * "f_alg.jl",         save_string(res_algFunc))
        write(dir * "J_y_alg.jl",       save_string(jac_algFunc))
        write(dir * "f_diff.jl",        save_string(res_diffFunc))

        @save dir * "J_sp.jl" J_y_sp θ_keys
    end

    if verbose println("Finished\n") end

    return Y0Func, res_algFunc, res_diffFunc, jacFunc, jac_algFunc, J_y_sp, θ_keys
end

function load_functions_forward_diff(p::AbstractModel)

    θ_sym, Y_sym, YP_sym, t_sym, SOC_sym, X_applied, γ_sym, p_sym, θ_keys = get_symbolic_vars(p)

    ## Y0 function
    Y0_sym = _symbolic_initial_guess(p_sym, SOC_sym, θ_sym, X_applied)

    ## batteryModel function
    res = _symbolic_residuals(p_sym, t_sym, Y_sym, YP_sym)

    θ_sym_slim, θ_keys_slim = get_only_θ_used_in_model(θ_sym, θ_keys, res, Y0_sym)

    Y0Func = build_function(Y0_sym, SOC_sym, θ_sym_slim, X_applied, fillzeros=false, checkbounds=false)[2]

    res_algFunc, res_diffFunc = _symbolic_initial_conditions_res(p, res, t_sym, Y_sym, YP_sym, θ_sym, θ_sym_slim)

    J_y_sp = _Jacobian_sparsity_pattern(p, res, Y_sym, YP_sym)[1]

    """
    Sometimes, the differential terms do not appear in the row corresponding to their value.
    this line ensures that there is still a spot for the ∂x/∂t term in the jacobian sparsity pattern/
    functon, otherwise there will be a mismatch in the sparse matrix
    """
    @inbounds for i in 1:p.N.diff
        if iszero(J_y_sp[i,i])
            J_y_sp[i,i] = 1e-100 # some arbitrarily small number so it doesn't get overwritten
        end
    end

    J_y_sp_alg = @inbounds J_y_sp[p.N.diff+1:end,p.N.diff+1:end]

    function build_color_jacobian_struct(J, f!, N, )
        colorvec = matrix_colors(J)

        Y_cache = zeros(Float64, N+1)

        _Y_cache  = zeros(Float64,p.N.tot)
        _YP_cache = zeros(Float64,p.N.tot)
        func = res_FD(f!, _Y_cache, _YP_cache, p.cache.θ_tot, p.N)

        jac_cache = ForwardColorJacCache(
            func,
            Y_cache,
            dx = zeros(N),
            colorvec = colorvec,
            sparsity = similar(J),
            )

        J_struct = jacobian_AD(func, J, jac_cache)
        return J_struct
    end

    initial_guess! = eval(Y0Func)
    f_alg!         = eval(res_algFunc)
    f_diff!        = eval(res_diffFunc)
    f! = function (res,t,Y,YP,θ_tot)
        f_diff!((@views @inbounds res[1:p.N.diff]),       t, Y, YP, θ_tot)
        f_alg!( (@views @inbounds res[p.N.diff+1:end-1]), t, Y, YP, θ_tot)
    end

    J_y!     = build_color_jacobian_struct(J_y_sp, f!, p.N.tot-1)
    J_y_alg! = build_color_jacobian_struct(J_y_sp_alg, f_alg!, p.N.alg-1)

    @assert size(J_y!.sp) == (p.N.tot-1,p.N.tot)
    @assert size(J_y_alg!.sp) == (p.N.alg-1,p.N.alg)

    return initial_guess!, f_alg!, f_diff!, J_y!, J_y_alg!, θ_keys_slim
end

function _symbolic_initial_guess(p::AbstractModel, SOC_sym, θ_sym, X_applied)

    Y0_sym = guess_init(p, X_applied)[1]

    deleteat!(Y0_sym, p.N.tot)

    return Y0_sym
end

function _symbolic_residuals(p::AbstractModel, t_sym, Y_sym, YP_sym)
    ## symbolic battery sol
    res = similar(Y_sym)
    residuals_PET!(res, t_sym, Y_sym, YP_sym, p)

    deleteat!(res, p.N.tot)

    return res
end

function _Jacobian_sparsity_pattern(p, res, Y_sym, YP_sym)
    sp_x = jacobian_sparsity(res, Y_sym)
    sp_xp = jacobian_sparsity(res, YP_sym)

    I,J,_ = findnz(sp_x .+ sp_xp)

    J_sp = sparse(I, J, 1.0, p.N.tot-1, p.N.tot)

    return J_sp, sp_x, sp_xp
end

function _symbolic_jacobian(p::AbstractModel=petlion(LCO);inds::T=1:p.N.tot) where T<:UnitRange{Int64}
    θ_sym, Y_sym, YP_sym, t_sym, SOC_sym, X_applied, γ_sym, p_sym, θ_keys = get_symbolic_vars(p)
    res = [_symbolic_residuals(p_sym, t_sym, Y_sym, YP_sym);Y_sym[end]]

    res    = @inbounds @views res[inds]
    Y_sym  = @inbounds @views Y_sym[inds]
    YP_sym = @inbounds @views YP_sym[inds]

    ## symbolic jacobian
    Jac_x  = sparsejacobian_multithread(res, Y_sym;  simplify=false)
    Jac_xp = sparsejacobian_multithread(res, YP_sym; simplify=false)

    Jac = Jac_x
    @inbounds for (i,j) in zip(Jac_xp.rowval,Jac_xp.colptr)
        Jac[i,j] += γ_sym*Jac_xp[i,j]
    end

    return Jac
end
function _symbolic_jacobian(p::AbstractModel, res, t_sym, Y_sym, YP_sym, γ_sym, θ_sym, θ_sym_slim; verbose=false)
    J_sp, sp_x, sp_xp = _Jacobian_sparsity_pattern(p, res, Y_sym, YP_sym)

    ## symbolic jacobian
    Jac_x  = sparsejacobian_multithread(res, Y_sym;  sp = sp_x,  simplify=false)
    Jac_xp = sparsejacobian_multithread(res, YP_sym; sp = sp_xp, simplify=false)

    Jac = Jac_x
    @inbounds for (i,j) in zip(Jac_xp.rowval,Jac_xp.colptr)
        Jac[i,j] += γ_sym*Jac_xp[i,j]
    end

    # building the sparse jacobian
    jacFunc = build_function(Jac.nzval, t_sym, Y_sym, YP_sym, γ_sym, θ_sym_slim)[2]

    return Jac, jacFunc, J_sp
end
function _symbolic_initial_conditions_res(p::AbstractModel, res, t_sym, Y_sym, YP_sym, θ_sym, θ_sym_slim)

    res_diff = res[1:p.N.diff]
    res_alg = res[p.N.diff+1:end]

    res_diffFunc = build_function(res_diff, t_sym, Y_sym, YP_sym, θ_sym_slim, fillzeros=false, checkbounds=false)[2]
    res_algFunc = build_function(res_alg,  t_sym, Y_sym, YP_sym, θ_sym_slim, fillzeros=false, checkbounds=false)[2]

    return res_algFunc, res_diffFunc
end
function _symbolic_initial_conditions_jac(p::AbstractModel, Jac, t_sym, Y_sym, YP_sym, γ_sym, θ_sym_slim)

    jac_alg = @inbounds Jac[p.N.diff+1:end,p.N.diff+1:end]

    jac_algFunc = build_function(jac_alg.nzval, t_sym, Y_sym, YP_sym, γ_sym, θ_sym_slim, fillzeros=false, checkbounds=false)[2]

    return jac_algFunc
end

function get_only_θ_used_in_model(θ_sym, θ_keys, X...)
    """
    Some of the parameters in θ may not be used in the sol. This returns
    a vector of the ones that are actually used and their sorted names
    """
    used_params_tot = eltype(θ_sym)[]
    @inbounds for _X in X, x in _X
        append!(used_params_tot, get_variables(x))
    end
    index_params = Int64[]

    dummy = BitArray{1}(undef, length(θ_sym))
    @inbounds for x in unique(used_params_tot)
        dummy .= isequal.(x, θ_sym)
        if sum(dummy) == 1
            index_param = findfirst(dummy)
            push!(index_params, index_param)
        end
    end
    sort!(index_params)

    θ_sym_slim = θ_sym[index_params]

    # for some reason, sometimes the above script does not remove all duplicates
    unique_ind = indexin(unique(θ_sym_slim), θ_sym_slim)
    θ_sym_slim = θ_sym_slim[unique_ind]

    index_params = index_params[unique_ind]

    θ_keys_slim = θ_keys[index_params]

    # Vectors of length n will have n-1 duplicates. This removes them
    unique_ind = indexin(unique(θ_keys_slim), θ_keys_slim)
    θ_keys_slim = θ_keys_slim[unique_ind]

    return θ_sym_slim, θ_keys_slim
end
@inline update_θ!(p::model) = update_θ!(p.cache.θ_tot,p.cache.θ_keys,p)
@inline function update_θ!(θ::Vector{Float64},keys::Vector{Symbol},p::model{jac}) where jac<:jacobian_symbolic
    θ_Dict = p.θ
    @inbounds for i in 1:length(θ)
        θ[i] = θ_Dict[keys[i]]
    end
    θ_Dict[:I1C] = calc_I1C(θ_Dict)
    return nothing
end
@inline function update_θ!(θ::Vector{Float64},keys::Vector{Symbol},p::model{jac}) where jac<:jacobian_AD
    θ_Dict = p.θ
    @inbounds for i in 1:length(θ)
        θ[i] = θ_Dict[keys[i]]
    end
    θ_Dict[:I1C] = calc_I1C(θ_Dict)

    @inbounds p.cache.θ_tot .= @views θ[1:length(p.cache.θ_tot)]
    return nothing
end

function get_symbolic_vars(p::AbstractModel;
    original_keys=nothing)
    if isnothing(original_keys)
        θ_keys = sort!(Symbol.(keys(p.θ)))
    else
        θ_keys = deepcopy(original_keys)
        all_keys = sort!(Symbol.(keys(p.θ)))

        @inbounds for key in all_keys
            if !(key ∈ θ_keys)
                push!(θ_keys, key)
            end
        end
    end

    @variables Y[1:p.N.tot], YP[1:p.N.tot], t, SOC, X_applied, γ
    @variables θ[1:length(θ_keys)]

    Y  = collect(Y)
    YP = collect(YP)
    θ  = collect(θ)

    p = model_skeleton([convert(_type,deepcopy(getproperty(p,field))) for (field,_type) in zip(fieldnames(model_skeleton),fieldtypes(model_skeleton))]...)
    p.opts.SOC = SOC

    @inbounds for (i,key) in enumerate(θ_keys)
        p.θ[key] = θ[i]
    end

    return θ, Y, YP, t, SOC, X_applied, γ, p, θ_keys
end

function sparsejacobian_multithread(ops::AbstractVector{<:Num}, vars::AbstractVector{<:Num};
    sp = jacobian_sparsity(ops, vars),
    simplify = true,
    multithread = true,
    )

    I,J,_ = findnz(sp)

    exprs = Vector{Num}(undef, length(I))

    if multithread
        @inbounds Threads.@threads for iter in 1:length(I)
            @inbounds exprs[iter] = expand_derivatives(Differential(vars[J[iter]])(ops[I[iter]]), simplify)
        end
    else
        @inbounds for iter in 1:length(I)
            @inbounds exprs[iter] = expand_derivatives(Differential(vars[J[iter]])(ops[I[iter]]), simplify)
        end
    end

    jac = sparse(I,J, exprs, length(ops), length(vars))
    return jac
end

remove_comments(str::String,x...) = remove_comments(collect(str),x...)
function remove_comments(str::Vector{Char},first::String="#=",last::String="=#")
    """
    An annoying aspect of Symbolics is that saving the output of `build_function` has tons of comments,
    severely adding to the memory size of the stored functions. This removes all comments
    """

    ind_first = find_in_string(str, first, 1)
    ind_last = find_in_string(str, last, length(last))

    @assert length(ind_first) == length(ind_last)

    ind_first = reverse(ind_first)
    ind_last  = reverse(ind_last)

    @inbounds for (l,u) in zip(ind_first,ind_last)
        index = l:u
        deleteat!(str,index)
    end
    return str
end

find_in_string(str::String,x...) = find_in_string(collect(str),x...)
function find_in_string(str::AbstractVector,x::String,I::T=1:length(x)) where T
    N = length(str)

    m = length(x)

    ind_vec = T[]
    @inbounds for i in 1:N - (m-1)
        index = i:(i+(m-1))
        y = @views @inbounds join(str[index])
        if y == x
            push!(ind_vec, (@inbounds index[I]))
        end
    end
    return ind_vec
end

function replace_repeated(str,x...)
    len_prev = length(str)
    len_new = 0
    @inbounds while len_prev != len_new
        len_prev = length(str)
        str = replace(str, x...)
        len_new = length(str)
    end
    return str
end

find_next(str::String,x...;kw...) = find_next(collect(str),x...;kw...)
function find_next(str,first,x::AbstractArray;itermax=900000)
    ind = first .+ (0:length(x)-1)
    iter = 0
    while (@views @inbounds str[ind]) != x
        ind = ind .+ 1
        iter += 1
        if ind[end] > length(str) error("Couldn't find the string") end
    end
    return ind
end
function find_next(str,first,x;itermax=900000)
    ind = first
    iter = 0
    while (@inbounds str[ind]) != x
        ind = ind .+ 1
        iter += 1
        if iter == itermax error(str[ind]) end
    end
    return ind
end

rearrange_if_statements(str::String) = rearrange_if_statements(collect(str))
function rearrange_if_statements(str::Vector{Char})
    inds = find_in_string(str, "if")
    @inbounds for ind_start in reverse(inds), _ in 1:4
        ind = find_next(str,ind_start[1] .+ (-1:0),['\n', ' '])
        while str[ind] == ['\n', ' ']
            deleteat!(str, ind[end])
        end

        str[ind[1]] = ';'
    end
    str = join(str)

    return str
end

convert_to_ifelse(str::String) = convert_to_ifelse(collect(str))
function convert_to_ifelse(str::Vector{Char})

    str = join(str)
    str = replace(str, "; else; " => ", ")
    str = replace(str, "; end" => ")")
    str = replace(str, ";" => ",")
    str = replace(str, "if" => "IfElse.ifelse(")

    return str
end


text:
function phi_galerk_modelA(phi, V, t, dt, lambda, r, g, k, xstep, T, D)
    noise = sqrt(T * dt / xstep) * sampler(size(phi))
    d = length(size(phi))
    R = CartesianIndices(phi)
    out = similar(phi)
    @. out = (phi - dt * (g * phi^2 + r + lambda * V[1]) * phi + noise)
    out = fft(out)
    for I in R
        out[I] = 1.0 / (1.0 + dt * sum([k[I[i]]^2 for i = 1:d])) * out[I]
    end
    out = ifft(out)
    @. phi = real(out)
end

#t =0.0
#r = 0.1
#g = 0.0
#D = 1.0
#T = 1.0
#lambda = 0.0
#@btime phi_galerk_modelA(phi, V, t, dt, lambda, r, g, k, xstep, T, D)

max_stars_repo_name: phyjonas/ImpuritySGPE

max_issues_repo_name: phyjonas/ImpuritySGPE

max_forks_repo_name: phyjonas/ImpuritySGPE


text:
# Command: Find Element from Element
# Refer to `src/01_structs/06_Element.jl`


text:
import ..UncertainValues.AbstractUncertainScalarKDE

"""
    resample(uv::UncertainScalarKDE)

Resample an uncertain value whose distribution is approximated using a
kernel density estimate once.
"""
resample(uv::AbstractUncertainScalarKDE) = rand(uv)

"""
    resample(uv::AbstractUncertainScalarKDE)

Resample an uncertain value whose distribution is approximated using a
kernel density estimate `n` times.
"""
resample(uv::AbstractUncertainScalarKDE, n::Int) = rand(uv, n)


text:
function maketablea2(dist::DataFrame;
    panela::Bool = true,
    panelb::Bool = true,
    )::Nothing

  #adjust the dataframe for the Harzmark Solomon stats
  hs::DataFrame = makehs(dist)

  #########################Panel a
  TA2_PANEL_A_ROWKEYS::Vector{Symbol} = [
    :ann,
    :interim,
    :ex,
    :annex,
    :ex40]

  TA2_PANEL_A_ROWLABEL::Dict = Dict(
    :ann=>"Announcement day",
    :interim=>"Interim",
    :ex=>"Ex-day",
    :annex=>"Ann-to-ex",
    :ex40=>"40 days after"
    )
  TA2_PANEL_A_COLNAMES::Vector{Vector{String}} = [["Return"]]

  #####################Panel b#######
  TA2_PANEL_B_ROWKEYS::Vector{Symbol} = [
    :constant,
    :yield,
    :r2,
    :N,]

  TA2_PANEL_B_ROWLABEL::Dict = Dict(
    :constant=>"Constant",
    :yield=>"Div yield",
    :r2=>"\$R^2\$",
    :N=>"N",
    )

  TA2_PANEL_B_COLKEYS::Vector{Symbol} = [
    :ann,
    :interim,
    :ex,
    :annex,
    :ex40
  ]

  TA2_PANEL_B_COLNAMES::Vector{Vector{String}} = [[
    "Announcement day", "Interim", "Ex-day", "Ann-to-ex", "40 days after"]]

  TA2_PANEL_B_COLFILTER = Dict(
    :ann=>:isann,
    :interim=>:isinterim,
    :ex=>:isex,
    :annex=>:isannex,
    :ex40=>:isex40
  )

  TA2_PANEL_B_NAME::String = "ta2panelb-$(REPLICATION_TYPE[])-$(DIST_TYPE[])-$(OUT_SUFFIX[])"

  TA2_PANEL_B_Σ = clusteredΣ!

  panelb && tablea2panelb(hs,
    rowkeys=TA2_PANEL_B_ROWKEYS,
    rowlabel=TA2_PANEL_B_ROWLABEL,
    tablename=TA2_PANEL_B_NAME,
    colkeys = TA2_PANEL_B_COLKEYS,
    colfilter = TA2_PANEL_B_COLFILTER,
    colnames = TA2_PANEL_B_COLNAMES,
    Σfunc = TA2_PANEL_B_Σ
  )
end

function tablea2panelb(hs;
  rowkeys::Vector{Symbol} = error("rowkeys is required"),
  rowlabel::Dict = error("rowlabel is required"),
  colkeys::Vector{Symbol} = error("colkeys is required"),
  colnames::Vector{Vector{String}} = error("colnames is required"),
  colfilter::Dict = error("colfilter is required"),
  tablename::String = error("table name is required"),
  Fret::Symbol = DIST_RET,
  Σfunc::Function = error("Σfunc is required"),
  tablepath::String = TABLE_PATH,
  decimals::Int = 3,)

  #the colmns vary by the views
  local specdfs::Vector{SubDataFrame} = Vector{SubDataFrame}()

  #all the regression sets are the same, we are just sub-sampling the data
  specs::FMSpecs = FMSpecs(length(colkeys))
  (i->push!(specs,
    yspec=Fret,
    xspec=:yield,
    xnames=[:intercept, :yield],
    clusterspec=[:permno, :day])).(1:length(colkeys))

  #println(describe(hs))

  specs.clusterspecs[3] = :permno
  #set the filters
  push!(specdfs, view(hs, hs[!, colfilter[colkeys[1]]], [Fret; :yield; :day; :permno; :eid]))
  push!(specdfs, view(hs, hs[!, colfilter[colkeys[2]]], [Fret; :yield; :day; :permno; :eid]))
  push!(specdfs, view(hs, hs[!, colfilter[colkeys[3]]], [Fret; :yield; :day; :permno; :eid]))
  push!(specdfs, view(hs, hs[!, colfilter[colkeys[4]]], [Fret; :yield; :day; :permno; :eid]))
  push!(specdfs, view(hs, hs[!, colfilter[colkeys[5]]], [Fret; :yield; :day; :permno; :eid]))

  #generate the regression content (will display in Finometrics template)
  contentrownames::Vector{String} = (s::Symbol->rowlabel[s]).(rowkeys[1:2])
  computeFMLMresults!(specdfs, specs, containsmissings=false)

  Ndescrows::Int = 2
  Ncols = length(colkeys)
  desccontent::Vector{Vector{String}} = (i->Vector{String}(undef, Ncols)).(1:Ndescrows)
  descrownames::Vector{String} = (s::Symbol->rowlabel[s]).(rowkeys[3:4])

  for c ∈ 1:Ncols
    desccontent[1][c] = "$(num2str(R²(specs.results[c]), decimals))"
    desccontent[2][c] = "$(length(unique(specdfs[c].eid)))"
  end

  #=for (i,lin) ∈ enumerate(specs.results)
    println("Iteration $i: ")
    display(Σfunc(lin))
  end=#

  tt::String = textable(specs.results, Σfunc, [:intercept, :yield],
    colnames=colnames,
    scaling=[1.,1.],
    descrownames=descrownames,
    desccontent=desccontent,
    decimaldigits=3)

  writetextable(tt, path=tablepath, outname="$tablename.tex")

  return nothing
end

#creates the neccessary fields to make the hs dataframe
function makehs(dist::DataFrame; Fret= DIST_RET)
  local hs::DataFrame

  hs=dist[dist.oneex .* dist.match,:]

#distcds::Vector{Vector{Int}} = (cd::Int->digits(cd)[4:-1:1]).(cumex.distcd) #(Reverse the order)
  hs = deleterows!(hs, (!).(completecases(view(hs, :, [:yield, Fret]))))

  #identify the announcement day
  N::Int = size(hs,1)
  hs.isann = Vector{Bool}(undef, N)
  hs.isinterim = Vector{Bool}(undef, N)
  hs.isex = Vector{Bool}(undef, N)
  hs.isannex = Vector{Bool}(undef, N)
  hs.isex40 = Vector{Bool}(undef, N)

  hs.tokeep = trues(N)
  @assert issorted(hs, [:eid, :day])

  hsss::GroupedDataFrame = groupby(hs, :eid)
  @mpar for i ∈ 1:length(hsss)
    shs::SubDataFrame = hsss[i]
    dvec::Vector{Int} = digits(shs[1,:distcd])[4:-1:1]
    idx::Int = searchsortedfirst(shs.day, -1)
    preexprice::Float64 = shs.day[idx] == -1 ? shs.price[idx] : -1.0

    #get the distribution codes and filter using the HS criteria
    if (dvec[1] ≠ 1) || (dvec[2] ≠ 2) || (dvec[3] > 5) || preexprice < 5.0
      shs.tokeep .= false
      continue
    end

    shs.isann .= shs.date .== shs.dclrdate
    shs.isex .= shs.day .== 0

    for r::DataFrameRow ∈ eachrow(shs)
      r.isinterim = (r.date > r.dclrdate) && (r.day < 0)
      r.isannex = r.isann || r.isinterim || r.isex
      r.isex40 = (r.day > 0) && (r.day ≤ 40)
    end
  end

  hs = hs[hs.tokeep,:]

  #println(size(hs))
  select!(hs, [:date, :isann, :isinterim, :isex, :isannex, :isex40,
    :permno, :day, Fret, :yield, :eid])

  #WARNING delete later
  #filter!(r::DataFrameRow->year(r.date) ≤ 2011, hs)
  #println(describe(hs))

  return hs
end


text:
__precompile__()

module SCS

if isfile(joinpath(dirname(@__FILE__), "..", "deps", "deps.jl"))
    include("../deps/deps.jl")
else
    error("SCS not properly installed. Please run Pkg.build(\"SCS\") and restart julia")
end

import Base.Libdl: RTLD_LAZY, RTLD_DEEPBIND, RTLD_GLOBAL, dlopen

function __init__()
    vnum = VersionNumber(SCS_version())
    # default binaries need access to Julia's lapack symbols
    if is_unix()
        dlopen(Base.liblapack_name, RTLD_LAZY|RTLD_DEEPBIND|RTLD_GLOBAL)
    end
    depsdir = realpath(joinpath(dirname(@__FILE__),"..","deps"))
    if !(vnum.major == 1 && vnum.minor in [1,2])
        error("Current SCS version installed is $(SCS_version()), but we require version 1.1.* or 1.2.*. On Linux and Windows, delete the contents of the `$depsdir` directory except for the files `build.jl` and `.gitignore`, then rerun Pkg.build(\"SCS\"). On OS X, run `using Homebrew; Homebrew.update()` in Julia.")
    end
    if is_linux() && vnum.major == 1 && vnum.minor == 1
        warn("Current SCS version installed is $(SCS_version()), but a newer version (1.2.*) is available. To upgrade, delete the contents of the `$depsdir` directory except for the files `build.jl` and `.gitignore`, then rerun Pkg.build(\"SCS\").")
    end
end

include("types.jl")
include("low_level_wrapper.jl")
include("high_level_wrapper.jl")
include("SCSSolverInterface.jl")  # MathProgBase interface

end # module


text:
abstract type Event end

#Captures an event and the time it takes place
struct TimedEvent
    event::Event
    time::Float64
end

#Comparison of two timed events - this will allow us to use them in a heap/priority-queue
isless(te1::TimedEvent,te2::TimedEvent) = te1.time < te2.time

#Generic events that we can always use
struct EndSimEvent <: Event end
struct LogStateEvent <: Event end

#External arrival to the first queue
struct ExternalArrivalEvent <: Event end

struct EndOfServiceAtQueueEvent <: Event
    q::Int #The index of the queue where service finished
end

struct InTransitEvent <: Event
    q::Int #The index of the destination queue
end


text:
# Preconditioners
# TODO explain a bit of theory, or refer to a source
# TODO homogenize with the terms interface

# A preconditioner is a type P that will get called as `Pk = P(basis, kpt)`
# and returns an object Pk that is expected to support
# ldiv! and mul! (both 2 and 3-arg versions)
# Additionally if the solver supports adaptive preconditioning
# it will call `precondprep!(P,X)` right before calling `ldiv!`

import LinearAlgebra: ldiv!
import LinearAlgebra: mul!

precondprep!(P, X) = P  # This API is also used in Optim.jl

"""
No preconditioning
"""
struct PreconditionerNone end
PreconditionerNone(basis, kpt) = I

"""
(simplified version of) Tetter-Payne-Allan preconditioning
↑ M.P. Teter, M.C. Payne and D.C. Allan, Phys. Rev. B 40, 12255 (1989).
"""
mutable struct PreconditionerTPA{T <: Real}
    basis::PlaneWaveBasis
    kpt::Kpoint
    kin::Vector{T}  # kinetic energy of every G
    mean_kin::Union{Nothing, Vector{T}}  # mean kinetic energy of every band
    default_shift::T # if mean_kin is not set by `precondprep!`, this will be used for the shift
end

function PreconditionerTPA(basis::PlaneWaveBasis{T}, kpt::Kpoint; default_shift=1) where T
    scaling = only([t for t in basis.model.term_types if t isa Kinetic]).scaling_factor
    kin = Vector{T}([scaling * sum(abs2, G + kpt.coordinate_cart)
                     for G in G_vectors_cart(kpt)] ./ 2)
    PreconditionerTPA{T}(basis, kpt, kin, nothing, default_shift)
end

@views function ldiv!(Y, P::PreconditionerTPA, R)
    if P.mean_kin === nothing
        ldiv!(Y, Diagonal(P.kin .+ P.default_shift), R)
    else
        Threads.@threads for n = 1:size(Y, 2)
            Y[:, n] .= P.mean_kin[n] ./ (P.mean_kin[n] .+ P.kin) .* R[:, n]
        end
    end
    Y
end
ldiv!(P::PreconditionerTPA, R) = ldiv!(R, P, R)
(Base.:\)(P::PreconditionerTPA, R) = ldiv!(P, copy(R))

# These are needed by eg direct minimization with CG
@views function mul!(Y, P::PreconditionerTPA, R)
    if P.mean_kin === nothing
        mul!(Y, Diagonal(P.kin .+ default_shift), R)
    else
        Threads.@threads for n = 1:size(Y, 2)
            Y[:, n] .= (P.mean_kin[n] .+ P.kin) ./ P.mean_kin[n] .* R[:, n]
        end
    end
    Y
end
(Base.:*)(P::PreconditionerTPA, R) = mul!(copy(R), P, R)

function precondprep!(P::PreconditionerTPA, X)
    P.mean_kin = [sum(real.(conj.(x) .* P.kin .* x)) for x in eachcol(X)]
end


text:
using Thebes, Luxor

function frame(scene, framenumber)
    background("lightblue1")
    setopacity(0.5)

    theta = rescale(framenumber, 1, scene.framerange.stop, 0, 2pi)
    perspective(500)

    eyepoint(Point3D(300cos(theta), 300sin(theta), 50 + 10sin(theta)))

    object  = make(Pyramid)
    scaleby!(object, 150, 150, 150)
    pin(object)
    axes3D(50)

end

width, height = (800, 800)
pyramidmovie = Movie(width, height, "pyramidmovie")

animate(pyramidmovie, [
    Scene(pyramidmovie, frame, 1:50)
    ],
    creategif=true,
    framerate=20,
    pathname="/tmp/pyramid-animation.gif")


text:
function particle_variable(data::PyArray{T,N}, units::Units{U,D}) where {T,N,D,U}
    Array(data*units)
end

function particle_variable(x::PyArray{T,N}, y::PyArray{T,N}, units::Units{U,D}) where {T,N,D,U}
    data = Point{2,T}.(x, y)
    Array(data*units)
end

function particle_variable(x::PyArray{T,N}, y::PyArray{T,N}, z::PyArray{T,N}, units::Units{U,D}) where {T,N,D,U}
    data = Point{3,T}.(x, y, z)
    Array(data*units)
end

function particle_variable(fr::PyObject, key::Symbol)
    py_obj = getproperty(fr, key)
    @debug "Reading data"
    py_data = py_obj.data
    if py_data isa Tuple
        @debug "Converting vector data"
        data = convert_it.(py_data)
        @debug "Reading units"
        units = get_units.(py_obj.units)
        @assert length(unique(units)) == 1
        particle_variable(data..., units[1])
    else
        @debug "Converting scalar data"
        data = convert_it(py_data)
        @debug "Reading units"
        units = get_units(py_obj.units)
        particle_variable(data, units)
    end
end


text:
using CompileMRI
using Test
using Pkg

try
    using RomeoApp
catch
    if !isdir(CompileMRI.pathof("RomeoApp"))
        CompileMRI.download_pkg("RomeoApp")
    else
        Pkg.develop("RomeoApp")
    end
    using RomeoApp
end

@testset "Compile ROMEO" begin
    include("romeo_test.jl")
    include("compile_romeo_test.jl")
end

try
    using ClearswiApp
catch
    if !isdir(CompileMRI.pathof("ClearswiApp"))
        CompileMRI.download_pkg("ClearswiApp")
    else
        Pkg.develop("ClearswiApp")
    end
    using ClearswiApp
end

@testset "Compile CLEARSWI" begin
    include("clearswi_test.jl")
    include("compile_clearswi_test.jl")
end


text:
using GaussianProcesses
using FileIO
using LinearAlgebra

# description = ""
# noise_cov_name = "noise_covariance_0001"
# file = "calibrate_convadj_to_lesbrary/loss_landscape_$(description).jld2"

# G = load(file)["G"]
# Φ1 = load(file)[noise_cov_name*"/Φ1"]
# Φ2 = load(file)[noise_cov_name*"/Φ2"]

# directory = "QuickCES/"
# isdir(directory) || mkdir(directory)

# pvalues = Dict(
#     :convective_κz => collect(0.075:0.025:1.025),
#     :background_κz => collect(0e-4:0.25e-4:10e-4),
# )

# ni = length(pvalues[:convective_κz])
# nj = length(pvalues[:background_κz])

Φ = Φ1 .+ Φ2

x = hcat([[pvalues[:convective_κz][i], pvalues[:background_κz][j]] for i = 1:ni, j = 1:nj]...)

not_nan_indices = findall(.!isnan.(Φ))
Φ = Φ[not_nan_indices]
x = x[:, not_nan_indices]

ces_directory = joinpath(directory, "QuickCES/")
isdir(ces_directory) || mkdir(ces_directory)

predict = trained_gp_predict_function(X, Φ)
using OceanLearning.Transformations: ZScore, normalize!
using Statistics

normalize!(Φ, ZScore(mean(Φ), std(Φ)))

ces_directory = joinpath(directory, "QuickCES/")
isdir(ces_directory) || mkdir(ces_directory)

# MZero
#  * Candidate solution
#     Final objective value:     -7.398974e+03

mZero = MeanZero()
# kern = Matern(5 / 2, [0.0 for _ in 1:ni*nj], 0.0) + SE(0.0, 0.0)
kern = Matern(5 / 2, [0.0, 0.0], 0.0)
gp = GP(x, Φ, mZero, kern, -2.0)

optimize!(gp)

xs = x[1, :]
ys = x[2, :]
Φ_predicted = predict(x)

using Plots
p = Plots.plot(gp)
Plots.savefig(p, joinpath(directory, "hello.png"))

plot_contour(eki, xs, ys, Φ_predicted, "GP_emulated", ces_directory; zlabel = "Φ", plot_minimizer=true, plot_scatters=false, title="GP-Emulated EKI Objective, Φ")
plot_contour(eki, xs, ys, Φ, "Original", ces_directory; zlabel = "Φ", plot_minimizer=true, plot_scatters=false, title="EKI Objective, Φ")

###
### Run MCMC on the emulated forward map
###

###
### Run MCMC on the true forward map
###


text:
# This file was generated, do not modify it. # hide
using Gnuplot
x, y = [0,1,2], [0,1,2]
Z = [10 10 10; 10 3 10; 10 2 10]
@gsp x y Z "w l lc 'red'"
save(term="png", output="./code/plt3d_ex1.png") # hide

max_stars_repo_name: lazarusA/gnuplot-examples

max_issues_repo_name: lazarusA/gnuplot-examples

max_forks_repo_name: lazarusA/gnuplot-examples


text:
module TestElkan

using ParallelKMeans
using Test
using Random

@testset "basic kmeans elkan" begin
    X = [1. 2. 4.;]
    res = kmeans(Elkan(), X, 1; n_threads = 1, tol = 1e-6, verbose = false)
    @test res.assignments == [1, 1, 1]
    @test res.centers[1] ≈ 2.3333333333333335
    @test res.totalcost ≈ 4.666666666666666
    @test res.converged

    res = kmeans(Elkan(), X, 2; n_threads = 1, init = [1.0 4.0], tol = 1e-6, verbose = false)
    @test res.assignments == [1, 1, 2]
    @test res.centers ≈ [1.5 4.0]
    @test res.totalcost ≈ 0.5
    @test res.converged
end

@testset "elkan no convergence yield last result" begin
    X = [1. 2. 4.;]
    res = kmeans(Elkan(), X, 2; n_threads = 1, init = [1.0 4.0], tol = 1e-6, max_iters = 1, verbose = false)
    @test !res.converged
    @test res.totalcost ≈ 0.5
end

@testset "elkan singlethread linear separation" begin
    Random.seed!(2020)

    X = rand(3, 100)
    res = kmeans(Elkan(), X, 3; n_threads = 1, tol = 1e-10, max_iters = 10, verbose = false)

    @test res.totalcost ≈ 14.16198704459199
    @test !res.converged
    @test res.iterations == 10
end

@testset "elkan multithread linear separation quasi two threads" begin
    Random.seed!(2020)

    X = rand(3, 100)
    res = kmeans(Elkan(), X, 3; n_threads = 2, tol = 1e-6, verbose = false)

    @test res.totalcost ≈ 14.16198704459199
    @test res.converged
end

end # module


text:
#=

  The BUGS Book
  Example 2.7.1 Surgery (continued): prediction

  Cf ~/jags/bugs_book_2_7_1.R
  """
  model {
    theta ~ dbeta(3, 27) # prior distribution
    Y ~ dbin(theta, 20) # sampling distribution
    P6 <- step(Y - 5.5) # =1 if y >= 6, 0 otherwise
  }

Summary Statistics
  parameters      mean       std   naive_se      mcse         ess      rhat
      Symbol   Float64   Float64    Float64   Float64     Float64   Float64

       theta    0.0987    0.0524     0.0003    0.0011   1799.2552    1.0049
           y    1.9620    1.6057     0.0080    0.0367   1271.8025    1.0081

Quantiles
  parameters      2.5%     25.0%     50.0%     75.0%     97.5%
      Symbol   Float64   Float64   Float64   Float64   Float64

       theta    0.0193    0.0588    0.0900    0.1299    0.2219
           y    0.0000    1.0000    2.0000    3.0000    6.0000

Distributions of variable (num:40)
0.00000 =>   38931  (0.973275)
1.00000 =>    1069  (0.026725)
  4.176741 seconds (17.21 M allocations: 923.018 MiB, 28.68% gc time)

=#

using Turing, StatsPlots, DataFrames
include("jl_utils.jl")

@model function bugs_book_2_7_1()

    theta ~ Beta(3,27) # prior distribution
    y ~ Binomial(20, theta) # sampling distribution
    p6 ~ Dirac(y >= 6)

end

model = bugs_book_2_7_1()

num_chns = 4

# chns = sample(model, Prior(), MCMCThreads(), 10_000, num_chns)
# chns = sample(model, MH(), 100_000)
# chns = sample(model, MH(
#                        # :alpha => Normal(2,sqrt(2)),
#                        # :beta => Normal(2,sqrt(2)),
#                        # :sigma => Gamma(2,2)
#                        ), MCMCThreads(), 40_000, num_chns)

# chns = sample(model, PG(15), MCMCThreads(), 1_000, num_chns)
# chns = sample(model, SMC(10_000), MCMCThreads(), 20_000, num_chns)
chns = sample(model, IS(), 10_000)

# Both HMC and NUTS give the following error:
# ERROR: LoadError: TaskFailedException:
# TypeError: in typeassert, expected Float64, got a value of type ForwardDiff.Dual{Nothing,Float64,3}
# chns = sample(model, HMC(0.1,5), MCMCThreads(), 10_000, num_chns) # Error
# chns = sample(model, NUTS(0.65), MCMCThreads(), 10_000, num_chns) # Error

display(chns)
# display(plot(chns))

show_var_dist_pct(chns,:p6 )
show_var_dist_pct(chns,:y )
# show_var_dist_pct(chns,:theta )

max_stars_repo_name: tias/hakank

max_issues_repo_name: tias/hakank

max_forks_repo_name: tias/hakank


text:
# 3d
function Base.log(R::RotationVec)
    x, y, z = params(R)
    return @SMatrix [0 -z  y
                     z  0 -x
                    -y  x  0]
end

function Base.log(R::Rotation{3})
    log(RotationVec(R))
end

# 2d
function Base.log(R::Angle2d)
    θ, = params(R)
    return @SMatrix [0 -θ
                     θ  0]
end

function Base.log(R::Rotation{2})
    log(Angle2d(R))
end


text:
function match!(
    crr::CorefResolver,
    params::MatchParams,
    )
    #TODO: erspace1 needs to be the smaller
    #after extensions and restarts this needs to be verified
    source_pair, target_pair = crr |> innerspaces

    crs = crr.crs
    crspace = crr.crspace

    n1 = source_pair.erspace1.height
    n2 = source_pair.erspace2.height
    w = source_pair.erspace1.width
    crspace.free_i[1:n1] .= true
    crspace.free_j[1:n2] .= true
    prepspace = crspace.prepspace
    prepw = size(prepspace)[2]

    _distfun(i::Int, j::Int) = dist(source_pair, crs, i, j)

    prepmatch!(_distfun, prepspace, n1, n2)

    match_count = 0
    best_match = Match()
    while match_count < params.needed_matches
        best_dist = Inf
        @inbounds for i in 1:n1
            !crspace.free_i[i] && continue
            new_match = prepspace[i, 1]
            new_match.dist > best_dist && continue

            if !crspace.free_j[new_match.ind]
                for k in 1:prepw-1
                    prepspace[i, k] = prepspace[i, k+1]
                end
                if prepspace[i, 1].ind == -1
                    (_j, _dist) = min_ind_of_free(x -> _distfun(i, x), n2, crspace.free_j)
                    prepspace[i, 1] = HalfMatch(_j, _dist)
                else
                    prepspace[i, prepw] = HalfMatch()
                end
                new_match = prepspace[i, 1]
            end

            if new_match.dist < best_dist
                best_dist = new_match.dist
                best_match = Match(i, new_match.ind, new_match.dist)
            end

        end
        crspace.free_i[best_match.i] = false
        crspace.free_j[best_match.j] = false
        match_count += 1
        for (_sarr, _tarr, ind) in zip(
            [source_pair.erspace1, source_pair.erspace2],
            [target_pair.erspace1, target_pair.erspace2],
            [best_match.i, best_match.j])

            for k in 1:w _tarr.space[match_count, k] = _sarr.space[ind, k] end
        end
    end
    for erspace in target_pair |> spaces
        erspace.height = match_count
        erspace.width = w
        erspace.entity_types = source_pair.erspace1.entity_types
    end
    target_pair.erspace1.height = match_count
    target_pair.erspace2.height = match_count
end

const MATCH_DIST = -10.0
const PATH_DISCOUNT = 0.8

function dist(crs::CorefSystem, name::String, i::Int, j::Int)
    if (i, j) in crs.dict_of_matches[name]  # TODO: from latent vars
        return MATCH_DIST
    end
    dfp = crs.esp_dict[name]
    normalizer = crs.normalizers[name]
    s = 0.0
    for (colname, mu, rho) in zip(names(dfp.df1), normalizer.means, normalizer.stds)
        s += @inbounds (col_dist(dfp.df1[i,colname], dfp.df2[j,colname]) - mu) / rho
    end
    s #/ length(normalizer.means)
end

function dist(space_pair::SpacePair, crs::CorefSystem, i::Int, j::Int)
    s = 0.0
    for l in 1:space_pair.erspace1.width
        name = space_pair.erspace1.entity_types[l]
        ent_ind_i = space_pair.erspace1.space[i, l]
        ent_ind_j = space_pair.erspace2.space[j, l]
        s += dist(crs, name, ent_ind_i, ent_ind_j) * PATH_DISCOUNT ^ (l-1)
    end
    s
end

function min_ind_of_free(fun::Function, n::Int, free_inds::Array{Bool})
    m = Inf
    j = -1
    for i in 1:n
        !free_inds[i] && continue
        d = fun(i)
        d < m && (j = i; m = d)
    end
    (j, m)
end

function prepmatch!(
    distfun::Function,
    prepspace::Array{HalfMatch, 2},
    n1::Int,
    n2::Int
    )
    prep_size = size(prepspace)[2]
    for i in 1:n1
        prepv = prepspace[i, :]
        prepmax = Inf
        fill!(prepv, HalfMatch())
        for j in 1:n2
            d = distfun(i, j)
            if d < prepmax
                new_match = HalfMatch(j, d)
                pind = searchsortedfirst(prepv, new_match)
                pind < prep_size && (prepv[pind+1:prep_size] = prepv[pind:prep_size-1])
                prepv[pind] = new_match
                prepmax = prepv[prep_size].dist
            end
        end
    end
end


text:
using EDFPlus
using Test

edfh = loadfile("EDFPlusTestFile.edf")

EDFPlus.version() != nothing
epoch_iterator(edfh, 1) != nothing
annotation_epoch_iterator(edfh, 1) != nothing
channeltimesegment(edfh, 1, 0.5, 1.5, true) != nothing
multichanneltimesegment(edfh, [1], 0.25, 1.25, false) != nothing
sg = physicalchanneldata(edfh, 2) != nothing
(fs = samplerate(edfh, 2)) != nothing


text:
include("new_io_tree.jl")
using X86Assembly

X86Assembly.representative(x::Integer) = x
function X86Assembly.representative(x::AbstractBitBundle)
  val::UInt8 = 0
  val
end
X86Assembly.representative(x::X86Assembly.Operand) =
  X86Assembly.Operand{UInt8, UInt32}(x.deref, X86Assembly.representative(x.base), X86Assembly.representative(x.scale),
                                              X86Assembly.representative(x.index), X86Assembly.representative(x.disp))
X86Assembly.representative(x::X86Assembly.AVX512Operand) =
  X86Assembly.AVX512Operand(
                X86Assembly.representative(x.mask),
                X86Assembly.representative(x.zero_combining),
                X86Assembly.representative(x.sae),
                X86Assembly.representative(x.er),
                X86Assembly.representative(x.rounding_mode))

const IntOrAbstract = Union{Integer, AbstractByte}
function X86Assembly.Operand(x::Bool, base::IntOrAbstract, scale::IntOrAbstract, index::IntOrAbstract, disp::T) where T
    X86Assembly.Operand{AbstractByte, T}(x, base, scale, index, disp)
end

function X86Assembly.readImmediate(tree::GenerativeIOTree, T)
  tree.hadImmediate = true
  tree.nbytes_read += sizeof(T)
  AbstractImmediate(sizeof(T))
end

using X86Assembly: compressed_prefix_status

include("tree_folding.jl")
include("codegen.jl")

function run!(projection)
  mode = X86Assembly.processor_modes(true, true)
  n = 0
  #p = Progress(1616850)
  # Create a compressed structure from two bytes. We'll have the decoder work
  # from that.
  root = IOTreeRoot()
  tree = GenerativeIOTree(root, AbstractTrees.firststate(Leaves(root)), 1, false, 0)
  while true
    a = UInt16(read(tree, UInt8)) << 8
    b = read(tree, UInt8)
    compressed_prefix = a | b
    if !X86Assembly.validate_compressed_prefix(compressed_prefix_status(0b0,compressed_prefix,0))
        next!(tree, projection(tree, nothing)) || break
        continue
    end
    first_non_prefix_byte = read(tree, UInt8)
    prefix = compressed_prefix_status(0b0, compressed_prefix, first_non_prefix_byte)
    data = X86Assembly._parseOpc(tree, prefix, mode)
    # Create a tree of non-prefix bytes
    if !next!(tree, projection(tree, data))
        break
    end
    n += 1
  end
  fold!(root)
  function byte_reader(level)
       level == 1 ? :((prefix.flags >> 8) % UInt8) :
       level == 2 ? :(prefix.flags % UInt8) :
       level == 3 ? :(prefix.first_non_prefix_byte) :
                    :(read(io, UInt8))
  end
  code = generate_decoder_for_tree(root, byte_reader)
  thunk = quote
    function (io)
        prefix = X86Assembly.fast_compressed_prefix_scanner(io)
        $code
    end
  end
  f = eval(thunk)
  @show f
  println("Number of instruction variants: $n")
  (tree, f)
end
non_prefix_bytes(tree, data) = data == nothing ? -1 : tree.nbytes_read - 3
(tree, f) = run!(non_prefix_bytes)

function debug_tree(tree, io)
    prefix = X86Assembly.fast_compressed_prefix_scanner(io)
    bytes = [(prefix.flags >> 8)%UInt8, prefix.flags % UInt8, prefix.first_non_prefix_byte]
    AbstractTrees._print_tree(STDOUT, tree, 20; withinds = true) do buf, node, inds
        an = tree
        for (level, ind) in enumerate(inds)
            an = c = children(an)[ind]
            if level > length(bytes)
                push!(bytes, read(io, UInt8))
            end
            if !(bytes[level] in c.mask)
                showcompact(buf, node)
                return
            end
        end
        Base.print_with_color(:yellow, sprint(showcompact,node))
    end
end

buf = IOBuffer([0xf3,0x0f,0x58,0xed])
debug_tree(tree.tree, buf)


text:
function compensated_summation(x::T, δ::T, ε::T) where {T}
    a = x
    ε = ε + δ
    x = a + ε
    ε = ε + (a - x)

    return (x, ε)
end

function compensated_summation(x::VT, δ::VT, ε::VT) where {T, VT <: AbstractVector{T}}
    a  = copy(x)
    ε  = ε .+ δ
    x̄  = a .+ ε
    ε̄ = ε .+ (a .- x̄)

    return (x̄, ε̄)
end

function compensated_summation!(x::VT, δ::VT, ε::VT) where {T, VT <: AbstractVector{T}}
    a  = copy(x)
    ε  = ε .+ δ
    x .= a .+ ε
    ε .= ε .+ (a .- x)

    return (x, ε)
end


text:
@testset "Simplex" begin
    s = Simplex(1, 2, 3)
    @test s[:index] == 0
    @test s[:values] == [1,2,3]
    s[:index] = 1
    @test s[:index] == 1

    s = Simplex([1,2,3])
    @test dim(s) == 2
    @test s[:values] == [1,2,3]
    @testset "Simplex faces" for (o,t) in zip(faces(s), [Simplex(2, 3), Simplex(1, 3), Simplex(1, 2)])
        @test o[:values] == t[:values]
    end

    ss = Set{Simplex}()
    push!(ss, s)
    push!(ss, s)
    @test length(ss) == 1

    function Base.isless(a::Vector{Int},b::Vector{Int})
        length(a) < length(b) && return true
        length(a) > length(b) && return false
        return any([aa < bb for (aa,bb) in zip(a,b)])
    end
    @test [1] < [1,3]
    @test [1,2,1] > [1,3]
    @test [1,2] < [1,3]
    s = Simplex([1,2], [1,3])
    @test dim(s) == 1

    X = hcat(zeros(3), Matrix(I,3,3))
    S = Simplex(1,2,3,4)
    @test volume(S,X) == 1/6
end

@testset "Cube" begin
    c = Cube{Vector{Int}}([0,0,0],Int[]) # 0-cube
    @test dim(c) == 0
    c = Cube{Vector{Int}}([0,0,0],Int[0,0,0]) # 0-cube
    @test dim(c) == 0
    c = Cube{Vector{Int}}([0,0,0],[1]) # 1-cube
    @test dim(c) == 1
    c = Cube{Vector{Int}}([0,0,0],[1,0,1]) # 2-cube
    @test dim(c) == 2
    c = Cube{Vector{Int}}([0,0,0],[1,1,1]) # 3-cube
    @test dim(c) == 3
    @test c[:values] == ([0,0,0],[1,1,1])

    c = Cube{Vector{Float64}}([0.,0.,0.],[.5,.0,.5]) # 2-cube
    @test volume(c) == 0.5^2
    c = Cube{Vector{Float64}}([0.,0.,0.],[.5,.5,.5]) # 2-cube
    @test volume(c) == 0.5^3
end


text:
### A Pluto.jl notebook ###
# v0.12.7

using Markdown
using InteractiveUtils

# ╔═╡ 6fd616f6-1ea1-11eb-3814-8bfb4a096c49
using Pkg, DrWatson

# ╔═╡ 6ffe2628-1ea1-11eb-24ea-57f985146a72
begin
	@quickactivate "StatisticsWithJuliaPlutoNotebooks"
	using Distributions, Random, StatsBase, DataFrames, Plots
	Random.seed!(1)
end;

# ╔═╡ ed174bc4-1ea0-11eb-1e2f-a32874cec549
md"## Listing 3.4"

# ╔═╡ 70179b62-1ea1-11eb-13ce-7fb4d54bad4b
begin
	dist = TriangularDist(4,6,5)
	N = 10^6
	data = rand(dist,N)
	yData=(data .- 5).^2
end

# ╔═╡ 0b18657a-1f96-11eb-0dbe-614426f96594
("Mean: $(mean(yData)), Variance: $(var(data))")

# ╔═╡ 0b18a3f0-1f96-11eb-1f2b-0f12d42055fa
begin
	p1 = histogram(data, xlabel="x", bins=80, normed=true, ylims=(0,1.1))
	p2 = histogram(yData, xlabel="y", bins=80, normed=true, ylims=(0,15))
	plot(p1, p2, ylabel="Proportion", size=(800, 400), legend=:none)
end

# ╔═╡ 70182cd8-1ea1-11eb-094d-8d0d49cf15f3
md"## End of listing 3.4"

# ╔═╡ Cell order:
# ╟─ed174bc4-1ea0-11eb-1e2f-a32874cec549
# ╠═6fd616f6-1ea1-11eb-3814-8bfb4a096c49
# ╠═6ffe2628-1ea1-11eb-24ea-57f985146a72
# ╠═70179b62-1ea1-11eb-13ce-7fb4d54bad4b
# ╠═0b18657a-1f96-11eb-0dbe-614426f96594
# ╠═0b18a3f0-1f96-11eb-1f2b-0f12d42055fa
# ╟─70182cd8-1ea1-11eb-094d-8d0d49cf15f3


text:
using Test
using ForestBiometrics

@test emc(75.0,80.0) == 14.023800000000001


text:
import Test
import PredictMD

y_true = [3, -0.5, 2, 7,]
y_pred = [2.5, 0.0, 2, 8,]
Test.@test(
    isapprox(
        PredictMD.r2_score(y_true, y_pred,),
        0.948;
        atol = 0.001,
        )
    )

y_true = [1,2,3,]
y_pred = [1,2,3,]
Test.@test(
    isapprox(
        PredictMD.r2_score(y_true, y_pred,),
        1.0;
        )
    )

y_true = [1,2,3,]
y_pred = [2,2,2,]
Test.@test(
    isapprox(
        PredictMD.r2_score(y_true, y_pred,),
        0.0;
        )
    )

y_true = [1,2,3,]
y_pred = [3,2,1,]
Test.@test(
    isapprox(
        PredictMD.r2_score(y_true, y_pred,),
        -3.0;
        )
    )


text:
# ------  Access

Base.show(io :: IO, e :: Endowment{T1}) where T1 =
    print(io, typeof(e),  " $(name(e)) with ",  marginal(e));

name(e :: Endowment{T1}) where T1  =  e.name;
label(e :: Endowment{T1}) where T1  = e.label;
marginal(e :: Endowment{T1}) where T1  =  e.marginal;
Base.eltype(e :: Endowment{T1}) where T1 = T1;

validate_draws(e :: Endowment{T1}, x :: Vector{T1}) where T1  =
    validate_draws(marginal(e), x);

# ------  For testing

function make_test_endowment_vector()
    return [Endowment(:normal, "Normal", NormalMarginal(2.0, 1.5)),
        Endowment(:lognormal, "LogNormal", LogNormalMarginal(1.0, 2.0, 0.8)),
        Endowment(:uniform, "Uniform", UniformMarginal(-1.5, 0.5)),
        Endowment(:percentile, "Percentile", PercentileMarginal{Float32}()),
        Endowment(:bounded,  "Bounded", BoundedMarginal(3.0, 4.0)),
        Endowment(:unbounded,  "Unbounded", UnboundedMarginal{Float16}()),
        Endowment(:beta, "Beta", BetaMarginal(-1.2, 3.4, 1.4, 3.2)),
        Endowment(:unknown, "Unknown", UnknownMarginal{Vector{Int}}())]
end

draw_test_endowments(e :: Endowment{T1}, n, rng :: AbstractRNG) where T1  =
    draw_test_endowments(marginal(e), n, rng);

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


text:
#=
epigraph of ∞-norm of real or complex vectors
(u, w) : u ≥ ‖Wᵢ‖, ∀i
no extended formulation needed
real case (polyhedral): u ≥ Wᵢ, u ≥ -Wᵢ
complex case: (u, re(Wᵢ), im(Wᵢ)) in EpiNormEucl

dual cone is epigraph of 1-norm
(u, w) : u ≥ ∑ᵢ ‖Wᵢ‖
extended formulation: ∃ λ, u ≥ Σᵢ λᵢ
real case (polyhedral): λᵢ ≥ Wᵢ, λᵢ ≥ -Wᵢ
complex case: (λᵢ, re(Wᵢ), im(Wᵢ)) in EpiNormEucl
=#

mutable struct EpiNormInf{D <: PrimDual, C <: RealCompF, E <: NatExt} <: Cache
    oa_s::Vector{AE}
    d::Int
    λ::Vector{VR}
    W_temp::Vector{C}
    function EpiNormInf{D, C, E}() where {D <: PrimDual, C <: RealCompF, E <: NatExt}
        return new{D, C, E}()
    end
end

function MOIPajarito.Cones.create_cache(
    oa_s::Vector{AE},
    cone::Hypatia.EpiNormInfCone{RealF, C},
    opt::Optimizer,
) where {C <: RealCompF}
    dim = MOI.dimension(cone)
    @assert dim == length(oa_s)
    d = (C == CompF ? div(dim - 1, 2) : dim - 1)
    E = (cone.use_dual ? nat_or_ext(opt, d) : Nat)
    D = primal_or_dual(cone.use_dual)
    cache = EpiNormInf{D, C, E}()
    cache.oa_s = oa_s
    cache.d = d
    cache.W_temp = zeros(C, d)
    return cache
end

# primal cone functions

function MOIPajarito.Cones.add_init_cuts(cache::EpiNormInf{Prim, RealF}, opt::Optimizer)
    # add variable bound
    u = cache.oa_s[1]
    JuMP.@constraint(opt.oa_model, u >= 0)
    opt.use_init_fixed_oa || return

    # add all cuts necessary to define the polyhedral cone (no OA needed)
    @views w = cache.oa_s[2:end]
    JuMP.@constraints(opt.oa_model, begin
        u .>= w
        u .>= -w
    end)
    return
end

function MOIPajarito.Cones.add_init_cuts(cache::EpiNormInf{Prim, CompF}, opt::Optimizer)
    # add variable bound
    u = cache.oa_s[1]
    JuMP.@constraint(opt.oa_model, u >= 0)
    opt.use_init_fixed_oa || return

    # add cuts (rt2, ±1, ±1) on (u, re(Wᵢ), im(Wᵢ))
    @views w = cache.oa_s[2:end]
    for i in 1:(cache.d)
        re_w = w[2i - 1]
        im_w = w[2i]
        JuMP.@constraints(opt.oa_model, begin
            rt2 * u - re_w - im_w >= 0
            rt2 * u - re_w + im_w >= 0
            rt2 * u + re_w - im_w >= 0
            rt2 * u + re_w + im_w >= 0
        end)
    end
    return
end

function MOIPajarito.Cones.get_subp_cuts(
    z::Vector{RealF},
    cache::EpiNormInf{Prim, <:RealCompF},
    opt::Optimizer,
)
    R = cache.W_temp
    @views vec_copyto!(R, z[2:end])
    cuts = AE[]
    for (i, R_i) in enumerate(R)
        abs(R_i) < 1e-9 && continue
        cut = _get_cut(R_i, i, cache, opt)
        push!(cuts, cut)
    end
    return cuts
end

function MOIPajarito.Cones.get_sep_cuts(
    s::Vector{RealF},
    cache::EpiNormInf{Prim, <:RealCompF},
    opt::Optimizer,
)
    # decomposed gradient cut on (u, Wᵢ) is (1, -Wsᵢ / ‖Wsᵢ‖)
    us = s[1]
    @views Ws = vec_copyto!(cache.W_temp, s[2:end])
    cuts = AE[]
    for (i, Ws_i) in enumerate(Ws)
        abs_i = abs(Ws_i)
        if abs_i < opt.tol_feas || us - abs_i > -opt.tol_feas
            continue
        end
        R_i = Ws_i / -abs_i
        cut = _get_cut(R_i, i, cache, opt)
        push!(cuts, cut)
    end
    return cuts
end

function _get_cut(R_i::RealF, i::Int, cache::EpiNormInf{Prim, RealF}, opt::Optimizer)
    # strengthened disaggregated cut on (u, Wᵢ) is (‖Rᵢ‖, Rᵢ)
    u = cache.oa_s[1]
    w_i = cache.oa_s[1 + i]
    return JuMP.@expression(opt.oa_model, abs(R_i) * u + R_i * w_i)
end

function _get_cut(R_i::CompF, i::Int, cache::EpiNormInf{Prim, CompF}, opt::Optimizer)
    # strengthened cut on (u, re(Wᵢ), im(Wᵢ)) is (‖Rᵢ‖, re(Rᵢ), im(Rᵢ))
    u = cache.oa_s[1]
    re_w = cache.oa_s[2i]
    im_w = cache.oa_s[2i + 1]
    (re_i, im_i) = reim(R_i)
    return JuMP.@expression(opt.oa_model, abs(R_i) * u + re_i * re_w + im_i * im_w)
end

# dual cone functions

function MOIPajarito.Cones.get_subp_cuts(
    z::Vector{RealF},
    cache::EpiNormInf{Dual},
    opt::Optimizer,
)
    R = cache.W_temp
    @views vec_copyto!(R, z[2:end])
    return _get_cuts(R, cache, opt)
end

function MOIPajarito.Cones.get_sep_cuts(
    s::Vector{RealF},
    cache::EpiNormInf{Dual},
    opt::Optimizer,
)
    us = s[1]
    Ws = cache.W_temp
    @views vec_copyto!(Ws, s[2:end])
    if us - norm(Ws, 1) > -opt.tol_feas
        return AE[]
    end

    # gradient cut is (1, (-Wsᵢ / ‖Wsᵢ‖)ᵢ)
    R = zero(Ws)
    for (i, Ws_i) in enumerate(Ws)
        if abs(Ws_i) > 1e-10
            R[i] = -Ws_i / abs(Ws_i)
        end
    end
    return _get_cuts(R, cache, opt)
end

# unextended formulation

function MOIPajarito.Cones.add_init_cuts(
    cache::EpiNormInf{Dual, RealF, Nat},
    opt::Optimizer,
)
    # add variable bound
    u = cache.oa_s[1]
    JuMP.@constraint(opt.oa_model, u >= 0)
    opt.use_init_fixed_oa || return

    # add cuts u ≥ ±e'w, u ≥ |wᵢ|
    @views w = cache.oa_s[2:end]
    JuMP.@constraints(opt.oa_model, begin
        u >= sum(w)
        u >= -sum(w)
        u .>= w
        u .>= -w
    end)
    return
end

function MOIPajarito.Cones.add_init_cuts(
    cache::EpiNormInf{Dual, CompF, Nat},
    opt::Optimizer,
)
    # add variable bound
    u = cache.oa_s[1]
    JuMP.@constraint(opt.oa_model, u >= 0)
    opt.use_init_fixed_oa || return

    # add cuts rt2 * u ≥ ±e'w and linearizations of u ≥ |wᵢ|
    @views w = cache.oa_s[2:end]
    JuMP.@constraints(opt.oa_model, begin
        rt2 * u >= sum(w)
        rt2 * u >= -sum(w)
    end)
    for i in 1:(cache.d)
        re_w = w[2i - 1]
        im_w = w[2i]
        JuMP.@constraints(opt.oa_model, begin
            rt2 * u - re_w - im_w >= 0
            rt2 * u - re_w + im_w >= 0
            rt2 * u + re_w - im_w >= 0
            rt2 * u + re_w + im_w >= 0
        end)
    end
    return
end

function _get_cuts(R::Vector{C}, cache::EpiNormInf{Dual, C, Nat}, opt::Optimizer) where {C}
    # strengthened cut on (u, W) is (‖R‖∞, R)
    # TODO extreme ray decomposition?
    p = norm(R, Inf)
    r = reinterpret(RealF, R)
    u = cache.oa_s[1]
    @views w = cache.oa_s[2:end]
    cut = JuMP.@expression(opt.oa_model, p * u + JuMP.dot(r, w))
    return [cut]
end

# extended formulation

function MOIPajarito.Cones.num_ext_variables(cache::EpiNormInf{Dual, <:RealCompF, Ext})
    return cache.d
end

function MOIPajarito.Cones.extend_start(
    cache::EpiNormInf{Dual, C, Ext},
    s_start::Vector{RealF},
    opt::Optimizer,
) where {C}
    if s_start[1] < 1e-9
        return zeros(cache.d)
    end
    w_start = s_start[2:end]
    λ_start = abs.(reinterpret(C, w_start))
    return λ_start
end

function MOIPajarito.Cones.setup_auxiliary(
    cache::EpiNormInf{Dual, <:RealCompF, Ext},
    opt::Optimizer,
)
    @assert cache.d >= 2
    u = cache.oa_s[1]
    λ = cache.λ = JuMP.@variable(opt.oa_model, [1:(cache.d)], lower_bound = 0)
    JuMP.@constraint(opt.oa_model, u >= sum(λ))
    return λ
end

function MOIPajarito.Cones.add_init_cuts(
    cache::EpiNormInf{Dual, RealF, Ext},
    opt::Optimizer,
)
    opt.use_init_fixed_oa || return

    # add all cuts necessary to define the polyhedral cone (no OA needed)
    @views w = cache.oa_s[2:end]
    λ = cache.λ
    JuMP.@constraints(opt.oa_model, begin
        λ .>= w
        λ .>= -w
    end)
    return
end

function MOIPajarito.Cones.add_init_cuts(
    cache::EpiNormInf{Dual, CompF, Ext},
    opt::Optimizer,
)
    opt.use_init_fixed_oa || return

    # add cuts (rt2, ±1, ±1) on (λᵢ, re(Wᵢ), im(Wᵢ))
    @views w = cache.oa_s[2:end]
    for (i, λ_i) in enumerate(cache.λ)
        re_w = w[2i - 1]
        im_w = w[2i]
        JuMP.@constraints(opt.oa_model, begin
            rt2 * λ_i - re_w - im_w >= 0
            rt2 * λ_i - re_w + im_w >= 0
            rt2 * λ_i + re_w - im_w >= 0
            rt2 * λ_i + re_w + im_w >= 0
        end)
    end
    return
end

function _get_cuts(
    R::Vector{<:RealCompF},
    cache::EpiNormInf{Dual, <:RealCompF, Ext},
    opt::Optimizer,
)
    cuts = AE[]
    for (i, R_i) in enumerate(R)
        abs(R_i) < 1e-9 && continue
        cut = _get_cut(R_i, i, cache, opt)
        push!(cuts, cut)
    end
    return cuts
end

function _get_cut(R_i::RealF, i::Int, cache::EpiNormInf{Dual, RealF, Ext}, opt::Optimizer)
    # strengthened disaggregated cut on (λᵢ, Wᵢ) is (‖Rᵢ‖, Rᵢ)
    λ_i = cache.λ[i]
    w_i = cache.oa_s[1 + i]
    return JuMP.@expression(opt.oa_model, abs(R_i) * λ_i + R_i * w_i)
end

function _get_cut(R_i::CompF, i::Int, cache::EpiNormInf{Dual, CompF, Ext}, opt::Optimizer)
    # strengthened disaggregated cut on (λᵢ, re(Wᵢ), im(Wᵢ)) is (‖Rᵢ‖, re(Rᵢ), im(Rᵢ))
    λ_i = cache.λ[i]
    re_w = cache.oa_s[2i]
    im_w = cache.oa_s[2i + 1]
    (re_i, im_i) = reim(R_i)
    return JuMP.@expression(opt.oa_model, abs(R_i) * λ_i + re_i * re_w + im_i * im_w)
end