text:
typealias Constraint Vector{Int}

"""
    computescaleandshiftfromboundingbox(boundingBox::Vector{Float64})

Computes the scale and horizontal/vertical shift from a given bounding box.
This establishes an aspect-ratio-conserving transformation from the given
bounding box into a square that fits comfortably into the outer bounding triangle
of the DelaunayMeshes. The square (blue) fills ten per cent of the maximally
inscribed square (dashed).

![inscribed](inscribedsquare.svg)

Specifically, if ``(c_x, c_y)`` denote the center of the bounding box and
``\Delta = \max((x_\mathrm{max} - x_\mathrm{min}), (x_\mathrm{max} - x_\mathrm{min}))``
is the maximum extent of the bounding box, then the coordinate transformation to
the new primed coordinates is given by
``(x', y') = 0.1 (x - c_x, y - c_y)/\Delta + (c_x', c_y')``,
where ``(c_x', c_y') = (1.5, 1.25)`` denotes the center of the inscribed square.
"""
function computescaleandshiftfromboundingbox(boundingBox::Vector{Float64})
  local scale = max(boundingBox[2]-boundingBox[1], boundingBox[4]-boundingBox[3])
  local cx = 0.5*(boundingBox[2]+boundingBox[1])
  local cy = 0.5*(boundingBox[4]+boundingBox[3])
  return scale, cx, cy
end

"""
Defines a mesh. A mesh consists of a bounding box providing bounds to the
vertex coordinates, constraints consisting of closed polygons and the
corresponding Delaunay tesselation.
"""
type Mesh
   boundingBox::Nullable{Vector{Float64}}
   scale::Nullable{Float64}
   shiftX::Nullable{Float64}
   shiftY::Nullable{Float64}

   tesselation::DelaunayMeshes.DelaunayTesselation
   constraints::Vector{Constraint}
   faceLocation::Vector{Bool}
   voronoiVertices::Nullable{Array{Float64, 2}}

   Mesh() = new(
     Nullable{Vector{Float64}}(),
     Nullable{Vector{Float64}}(),
     Nullable{Vector{Float64}}(),
     Nullable{Vector{Float64}}(),
     DelaunayMeshes.DelaunayTesselation(),
     Vector{Constraint}(),
     [true, false],
     Nullable{Array{Float64, 2}}()
   )
end

"""
    getvoronoivertices(mesh::Mesh)

Computes the Voronoi vertices of the given DelaunayMeshes.

# Remarks
* The results are cached but need to be re-computed whenever a new vertex or
  a new constraint is pushed into the DelaunayMeshes.
"""
function getvoronoivertices(mesh::Mesh)
  if (!isnull(mesh.voronoiVertices))
    return mesh.voronoiVertices
  end

  mesh.voronoiVertices = fill(Inf, size(mesh.tesselation.faces, 1), 2)

  # go through all dual edges and compute barycenter of corresponding face
  for ei in 2:2:length(mesh.tesselation.edges)
    local curFace = mesh.tesselation.edges[ei].origin
    if get(mesh.voronoiVertices)[curFace, :] == [Inf, Inf]
      local curFaceEdge = rot(ei)
      local curFaceVertexIndices =
       lface_vertices(mesh.tesselation, curFaceEdge)
      local curFaceVertexCoordinates =
       map((i::Int) -> mesh.tesselation.vertices[i], curFaceVertexIndices)
      local baryCenter = sum(curFaceVertexCoordinates)/3.
      get(mesh.voronoiVertices)[curFace, :] = baryCenter
    end
  end

  return get(mesh.voronoiVertices)
end

function setboundingbox(mesh::Mesh, boundingBox::Vector{Float64})
  @assert isnull(mesh.boundingBox)
  mesh.boundingBox = boundingBox
  mesh.scale, mesh.shiftX, mesh.shiftY =
   computescaleandshiftfromboundingbox(get(mesh.boundingBox))
end

function scalepoints(mesh::Mesh, points::Array{Float64, 2})
  local scaledX = 0.05*(points[:, 1] - get(mesh.shiftX))/get(mesh.scale) + 1.5
  local scaledY = 0.05*(points[:, 2] - get(mesh.shiftY))/get(mesh.scale) + 1.25
  local scaledPoints = [scaledX scaledY]

  return scaledPoints
end

function unscalepoints(mesh::Mesh, points::Array{Float64, 2})
  local unScaledX = (points[:, 1]  - 1.5) * get(mesh.scale)/0.05 + get(mesh.shiftX)
  local unScaledY = (points[:, 2] - 1.25) * get(mesh.scale)/0.05 + get(mesh.shiftY)
  local unScaledPoints = [unScaledX unScaledY]

  return unScaledPoints
end

function push_scaled!(mesh::Mesh, scaledPoints::Array{Float64, 2})
  push!(mesh.tesselation, scaledPoints)

  # mark faces as interior
  local noMissingFaceLocations =
   length(mesh.tesselation.faces) - length(mesh.faceLocation)
  local truth = fill(true, noMissingFaceLocations)
  mesh.faceLocation = [mesh.faceLocation; truth ]

  # invalidate voronoi vertex cache
  mesh.voronoiVertices = Nullable{Array{Float64, 2}}()
end

"""
    push!(mesh::Mesh, points::Array{Float64, 2})

Inserts a number of vertices into the mesh.

# Remarks
* If the mesh was initialized without a bounding box, it is determined by
  the first set of vertices that are pushed into the triangualation. All
  consecutive push operations cannot insert vertices that exceed the bounding
  box limits
* New faces that are created by a push operation are marked as interior by
  default, i.e. it is forbidden to push vertices into an exterior region.
"""
function push!(mesh::Mesh, points::Array{Float64, 2})
  local minX = minimum(points[:,1])
  local maxX = maximum(points[:,1])
  local minY = minimum(points[:,2])
  local maxY = maximum(points[:,2])

  if isnull(mesh.boundingBox)
    setboundingbox(mesh, [minX, maxX, minY, maxY])
  else
    @assert (minX >= get(mesh.boundingBox)[1])
    @assert (maxX <= get(mesh.boundingBox)[2])
    @assert (minY >= get(mesh.boundingBox)[3])
    @assert (maxY <= get(mesh.boundingBox)[4])
  end

  local scaledPoints = scalepoints(mesh, points)
  push_scaled!(mesh, scaledPoints)
end

"""
    getfacesinsideregion(mesh::Mesh, ei::Int)

Finds all faces that are located inside the same region as the face identified
by `org(ei)`.

# Remarks
* The algorithm performs a breadth-first search over all adjacent triangles
  while omitting triangles that requires crossing a boundary.
"""
function getfacesinsideregion(mesh::Mesh, ei::Int)
  local foundFaces = Vector{Int}()
  local queuedEdges = [ei]
  while (length(queuedEdges) > 0)
    local currentEdge = shift!(queuedEdges)
    local currentFace = org(mesh.tesselation, currentEdge)

    # do nothing if it was already discovered
    if (currentFace ∉ foundFaces)
      push!(foundFaces, currentFace)

      # insert adjacent faces if they do not require boundary crossing
      local ec = currentEdge
      if (!isconstraint(mesh.tesselation, rot(ec)))
        push!(queuedEdges, sym(ec))
      end
      ec = onext(mesh.tesselation, ec)
      if (!isconstraint(mesh.tesselation, rot(ec)))
        push!(queuedEdges, sym(ec))
      end
      ec = onext(mesh.tesselation, ec)
      if (!isconstraint(mesh.tesselation, rot(ec)))
        push!(queuedEdges, sym(ec))
      end
    end
  end

  return foundFaces
end

"""
    addconstraint!(mesh::Mesh, vertices::Vector{Int})

Adds a constraint, given as a list of existing vertices, into the DelaunayMeshes.

# Remarks
* Constraints are closed polygons. A constraint edge from the last vertex to
  the first vertex is automatically inserted to ensure closure.
* All faces interior of the constrained region are marked as "exterior" to the
  DelaunayMeshes. A face is interior to the constrained region if it is located
  to the right of the given constraint polygon (i.e. the constraint vertices
  circle the interior region counter-clockwise).
* Constraint edges must not intersect. For performance reasons, this is not
  checked.
* Any vertex of the triangulation can intersect either zero or two constraint edges.
"""
function addconstraint!(mesh::Mesh, vertices::Vector{Int})
  # insert constraints between vertices and last one to close
  for i=1:length(vertices)-1
    insertconstraint!(mesh.tesselation, vertices[i], vertices[i+1])
  end
  local ci = insertconstraint!(mesh.tesselation, vertices[end], vertices[1])

  # find vertices to the right of the constraint
  local facesInside = getfacesinsideregion(mesh, rot(ci))
  mesh.faceLocation[facesInside] = false

  # invalidate voronoi vertex cache
  mesh.voronoiVertices = Nullable{Array{Float64, 2}}()
end

function istrianglevalid(mesh::Mesh, x::Tuple{Int, Int, Int})
  local rei = rot(x[1])
  local faceI = dest(mesh.tesselation, rei)
  return mesh.faceLocation[faceI]
end

function verticesfromtriangle(mesh::Mesh, x::Tuple{Int, Int, Int})
  local v1 = org(mesh.tesselation, x[1])
  local v2 = org(mesh.tesselation, x[2])
  local v3 = org(mesh.tesselation, x[3])

  return [v1 v2 v3]
end

"""
    convert(::Type{DiffEqPDEBase.SimpleFEMMesh}, mesh::Mesh)

Converts the given mesh to a format that can be used by the FEM-solvers.
"""
function convert(::Type{DiffEqPDEBase.SimpleFEMMesh}, mesh::Mesh)
  local allTriangles = DelaunayMeshes.getalltriangles(mesh.tesselation)
  local validTriangles = filter((x)->istrianglevalid(mesh,x), allTriangles)
  local triangles = map((x)->verticesfromtriangle(mesh,x), validTriangles)

  local redTriangles = reduce(vcat, triangles)
  local redVertices = reduce(hcat, mesh.tesselation.vertices)'

  return DiffEqPDEBase.SimpleFEMMesh(mesh.tesselation.vertices, redTriangles)
end


text:
using Documenter, Literate, UnitfulRecipes

# generate tutorials and how-to guides using Literate
src = joinpath(@__DIR__, "src")
lit = joinpath(@__DIR__, "lit")
notebooks = joinpath(src, "notebooks")

ENV["GKS_ENCODING"] = "utf-8"

execute = true # Set to true for executing notebooks and documenter!
nb = true      # Set to true to generate the notebooks
for (root, _, files) in walkdir(lit), file in files
    splitext(file)[2] == ".jl" || continue
    ipath = joinpath(root, file)
    opath = splitdir(replace(ipath, lit=>src))[1]
    Literate.markdown(ipath, opath, documenter = execute)
    nb && Literate.notebook(ipath, notebooks, execute = execute)
end

# Documentation structure
ismd(f) = splitext(f)[2] == ".md"
pages(folder) = [joinpath(folder, f) for f in readdir(joinpath(src, folder)) if ismd(f)]

makedocs(
    sitename="UnitfulRecipes.jl",
    doctest = false, # TODO guessing I should remove that when actually deploying?
    # options
    modules = [UnitfulRecipes],
    # organisation
    pages = Any[
        "Home" => "index.md",
        "Examples" => pages("examples")
    ]
)

# Deploy
deploydocs(
    repo = "github.com/jw3126/UnitfulRecipes.jl.git",
    push_preview = true
)

#=
To edit locally, make sure execute is set to false and run

using LiveServer
servedocs(literate=joinpath("docs", "lit"), doc_env=true)

from the root of the package in development
=#

ext: jl


text:
# using BioMedQuery.Processes
# using BioMedQuery.Entrez
# using BioMedQuery.UMLS
# using JLD
# using MySQL
# using BCBIVizUtils
# using Colors
using JLD
using PubMedMiner

"""
    epilepsy_comorbidities(;overwrite = true)

PuMed/Medline comorbidities for Epilepsy

Code by: Isabel Restrepo
BCBI - Brown University
Version: Julia 0.6

Before running, configure your MySQL onnection by adding to ~/.juliarc.jl
the following variables:

@everywhere ENV["PUBMEDMINER_DB_HOST"]="db_host"
@everywhere ENV["PUBMEDMINER_DB_USER"]="your_user"
@everywhere ENV["PUBMEDMINER_DB_PSSWD"]="your_password"

You alseo need to configure the email addess associated with you NCBI account,
and USER and PASSWD associated with your UMLS/UTS account.
"""
function epilepsy_semantic_occurrences(;overwrite=true)

    # Settings
    mh = "Epilepsy"
    concept = "Disease or Syndrome"

    info("----------------Start: umls_semantic_occurrences")
    save_semantic_occurrences(mh, concept; overwrite = overwrite)
end

# julia> @time epilepsy_semantic_occurrences(overwrite=true)
# INFO: ----------------Start: umls_semantic_occurrences
# INFO: 66720 Articles related to MH:Epilepsy
# INFO: ----------------------------------------
# INFO: Start all articles
# INFO: Using concept table: MESH_T047
# 211.869394 seconds (8.32 M allocations: 281.440 MiB, 0.04% gc time)

"""
epilepsy_comorbidities(;overwrite = true)

PuMed/Medline comorbidities for Epilepsy

Code by: Isabel Restrepo
BCBI - Brown University
Version: Julia 0.6

Before running, configure your MySQL onnection by adding to ~/.juliarc.jl
the following variables:

@everywhere ENV["PUBMEDMINER_DB_HOST"]="db_host"
@everywhere ENV["PUBMEDMINER_DB_USER"]="your_user"
@everywhere ENV["PUBMEDMINER_DB_PSSWD"]="your_password"

You alseo need to configure the email addess associated with you NCBI account,
and USER and PASSWD associated with your UMLS/UTS account.
"""
function colonic_neoplasms_semantic_occurrences(;overwrite=true)

    # Settings
    mh = "Colonic Neoplasms"
    concepts = ("Disease or Syndrome", "Mental or Behavioral Dysfunction", "Neoplastic Process");
    info("----------------Start: umls_semantic_occurrences")
    save_semantic_occurrences(mh, concepts...; overwrite = overwrite)
end

# @time colonic_neoplasms_semantic_occurrences(;overwrite=true)
# julia -p 2 comorbitidies_local_medline.jl
# INFO: ----------------Start: umls_semantic_occurrences
# INFO: 63304 Articles related to MH:Colonic Neoplasms
# INFO: ----------------------------------------
# INFO: Start all articles
# INFO: Using concept table: MESH_T047
# Disease or Syndrome: 271.139750 seconds (7.96 M allocations: 271.182 MiB, 0.03% gc time)
# Mental or Behavioral Dysfunction: Short...
# Neoplastic Process: 1120.822182 seconds (8.06 M allocations: 275.496 MiB, 0.01% gc time)


text:
using Sixel
using Test
using ImageCore, IndirectArrays, OffsetArrays
using ImageQualityIndexes
using LinearAlgebra
using FileIO, TestImages

sixel_output = Sixel.is_sixel_supported()
sixel_output || @info "Current terminal does not support sixel format sequence. Display tests to stdout will be marked as broken."
function test_sixel_display(f)
    if sixel_output
        @test_nowarn f()
    else
        @test_broken false
    end
end

@testset "Sixel.jl" begin
    include("backend/libsixel.jl")
    # include("fileio.jl") # need upstream registration
end


text:
data = readline()

function solve_1()
    count(==('('), data) - count(==(')'), data)
end

function solve_2()
    pos = 0
    for (i, c) in enumerate(data)
        pos += c == '(' ? 1 : -1
        if pos == -1
            return i
        end
    end
end

solve_1() |> println
solve_2() |> println

max_stars_repo_name: pLOPeGG/JuliaOfCode

max_issues_repo_name: pLOPeGG/JuliaOfCode

max_forks_repo_name: pLOPeGG/JuliaOfCode


text:
using StanModels, CSV

df = CSV.read(joinpath(@__DIR__, "..", "..", "data", "chimpanzees.csv"), delim=';');

# Define the Stan language model

m10_2s = "
data{
    int N;
    int pulled_left[N];
    int prosoc_left[N];
}
parameters{
    real a;
    real bp;
}
model{
    vector[N] p;
    bp ~ normal( 0 , 10 );
    a ~ normal( 0 , 10 );
    for ( i in 1:N ) {
        p[i] = a + bp * prosoc_left[i];
        p[i] = inv_logit(p[i]);
    }
    pulled_left ~ binomial( 1 , p );
}
";

# Define the Stanmodel and set the output format to :mcmcchains.

sm = SampleModel("m10.2s", m10_2s);

# Input data for cmdstan

m10_2_data = Dict("N" => size(df, 1),
"pulled_left" => df[!, :pulled_left], "prosoc_left" => df[!, :prosoc_left]);

# Sample using cmdstan

(sample_file, log_file)= stan_sample(sm, data=m10_2_data);

# Result rethinking

rethinking = "
   mean   sd  5.5% 94.5% n_eff Rhat
a  0.04 0.12 -0.16  0.21   180 1.00
bp 0.57 0.19  0.30  0.87   183 1.01
";

# Describe the draws

if !(sample_file == nothing)
  chn = read_samples(sm)
  describe(chn)
end


text:
using StaticArrays, EllZeroTrendFiltering
using Test

#SNP example
N = 400
data = snp500_data()[1:N]

M = 10
@testset "Examples lazy=$lazy" for lazy in [true, false]
    Ivec, Yvec, fvec = fit_pwl_constrained(data, M, lazy=lazy)

    #Make sure nothing crashes
    @test Ivec[5] == [1,78,146,162,370,400]
    @test Yvec[5] ≈ [7.180172223381005, 7.215552904744207, 7.160379047774212,
                     7.2454671184257435, 7.302625286124225, 7.215712991948725]
    @test fvec[5] ≈ 0.2057753829867579

    @test Ivec[M] == [1,64,68,144,164,204,241,252,281,366,400]
    @test Yvec[M] ≈ [7.192869876795329,7.187186703983764,7.234241678532294,7.156765323782069,
                  7.249687632432071,7.274247888087285,7.211190688270137,7.324697309669883,
                  7.258287507540115,7.316043144146352,7.215020359632374]
    @test fvec[M] ≈ 0.11431141407229006

    #Test regularize discrete
    # should generate solution with 10 segments
    I, Y, f = fit_pwl_regularized(data, 0.02, lazy=lazy)
    @test length(I) == 9
    @test I == Ivec[8]
    @test Y ≈ Yvec[8]     rtol=sqrt(eps())
    @test f ≈ fvec[8]     rtol=sqrt(eps())

    #Sin example
    g_(x) = sin(x) + 0.5sin(3.5x) + 0.5sin(5.1x)
    local t # TODO, can remove in julia 1.0?
    t = range(0, stop=2π, length=201)

    ζ = 0.1
    I, Y, cost = fit_pwl_regularized(g_, t, ζ, lazy=lazy)

    @test I == [1, 87, 107, 129, 147, 201]
    @test cost ≈ 0.35055983342102515

    # Compare to constrained
    Ivec, Yvec, fvec = fit_pwl_constrained(g_, t, 5, lazy=lazy)
    @test I == Ivec[5]
    @test Y ≈ Yvec[5]       rtol=sqrt(eps())
    @test cost ≈ fvec[5]       rtol=sqrt(eps())

    ζ = 0.002
    I, Y, cost = fit_pwl_regularized(g_, t, ζ, lazy=lazy)

    @test I == [1, 9, 15, 33, 52, 68, 87, 104, 111, 125, 132, 146, 153, 170, 188, 201]
    @test cost ≈ 0.006061712727833957
end


text:
@doc doc"""
    usinc(θ::Real)

Unnormalized version of `sinc` function, i.e.
$\operatorname{usinc}(\theta) = \frac{\sin(\theta)}{\theta}$. This is
equivalent to `sinc(θ/π)`.
"""
@inline usinc(θ::Real) = θ == 0 ? one(θ) : isinf(θ) ? zero(θ) : sin(θ) / θ

@doc doc"""
    usinc_from_cos(x::Real)

Unnormalized version of `sinc` function, i.e.
$\operatorname{usinc}(\theta) = \frac{\sin(\theta)}{\theta}$, computed from
$x = cos(\theta)$.
"""
@inline function usinc_from_cos(x::Real)
    if x >= 1
        return one(x)
    elseif x <= -1
        return zero(x)
    else
        return sqrt(1 - x^2) / acos(x)
    end
end

"""
    eigen_safe(x)

Compute the eigendecomposition of `x`. If `x` is a `StaticMatrix`, it is
converted to a `Matrix` before the decomposition.
"""
@inline eigen_safe(x; kwargs...) = eigen(x; kwargs...)

@inline function eigen_safe(x::StaticMatrix; kwargs...)
    s = size(x)
    E = eigen!(Matrix(parent(x)); kwargs...)
    return Eigen(SizedVector{s[1]}(E.values), SizedMatrix{s...}(E.vectors))
end

"""
    log_safe(x)

Compute the matrix logarithm of `x`. If `x` is a `StaticMatrix`, it is
converted to a `Matrix` before computing the log.
"""
@inline log_safe(x) = log(x)

@inline function log_safe(x::StaticMatrix)
    s = Size(x)
    return SizedMatrix{s[1],s[2]}(log(Matrix(parent(x))))
end

"""
    select_from_tuple(t::NTuple{N, Any}, positions::Val{P})

Selects elements of tuple `t` at positions specified by the second argument.
For example `select_from_tuple(("a", "b", "c"), Val((3, 1, 1)))` returns
`("c", "a", "a")`.
"""
@generated function select_from_tuple(t::NTuple{N, Any}, positions::Val{P}) where {N, P}
    for k in P
        (k < 0 || k > N) && error("positions must be between 1 and $N")
    end
    return Expr(:tuple, [Expr(:ref, :t, k) for k in P]...)
end

"""
    ziptuples(a, b)

Zips tuples `a` and `b` in a fast, type-stable way. If they have different
lengths, the result is trimmed to the length of the shorter tuple.
"""
@generated function ziptuples(a::NTuple{N,Any}, b::NTuple{M,Any}) where {N,M}
    ex = Expr(:tuple)
    for i = 1:min(N, M)
        push!(ex.args, :((a[$i], b[$i])))
    end
    ex
end

"""
    ziptuples(a, b, c)

Zips tuples `a`, `b` and `c` in a fast, type-stable way. If they have different
lengths, the result is trimmed to the length of the shorter tuple.
"""
@generated function ziptuples(a::NTuple{N,Any}, b::NTuple{M,Any}, c::NTuple{L,Any}) where {N,M,L}
    ex = Expr(:tuple)
    for i = 1:min(N, M, L)
        push!(ex.args, :((a[$i], b[$i], c[$i])))
    end
    ex
end

"""
    ziptuples(a, b, c, d)

Zips tuples `a`, `b`, `c` and `d` in a fast, type-stable way. If they have
different lengths, the result is trimmed to the length of the shorter tuple.
"""
@generated function ziptuples(a::NTuple{N,Any}, b::NTuple{M,Any}, c::NTuple{L,Any}, d::NTuple{K,Any}) where {N,M,L,K}
    ex = Expr(:tuple)
    for i = 1:min(N, M, L, K)
        push!(ex.args, :((a[$i], b[$i], c[$i], d[$i])))
    end
    ex
end

"""
    ziptuples(a, b, c, d, e)

Zips tuples `a`, `b`, `c`, `d` and `e` in a fast, type-stable way. If they have
different lengths, the result is trimmed to the length of the shorter tuple.
"""
@generated function ziptuples(a::NTuple{N,Any}, b::NTuple{M,Any}, c::NTuple{L,Any}, d::NTuple{K,Any}, e::NTuple{J,Any}) where {N,M,L,K,J}
    ex = Expr(:tuple)
    for i = 1:min(N, M, L, K, J)
        push!(ex.args, :((a[$i], b[$i], c[$i], d[$i], e[$i])))
    end
    ex
end

# TODO: make a better implementation for StaticArrays
function LinearAlgebra.eigvals(A::StaticArray, B::StaticArray; kwargs...)
    return eigvals(Array(A), Array(B); kwargs...)
end


text:
using Distributions: Normal, logpdf
using LinearAlgebra: ldiv!, diag, inv, tr

########################
# Leave-one-out        #
########################

function predict_LOO(Σ::AbstractPDMat, alpha::AbstractVector{<:Real}, y::AbstractVector{<:Real})
    invΣ = inv(Σ)
    σi2 = 1 ./ diag(invΣ)
    μi = -alpha .* σi2 .+ y
    return μi, σi2
end
"""
    predict_LOO(gp::GPE)

Leave-one-out cross-validated predictions.
Returns predictions of yᵢ given all other observations y₋ᵢ,
as a vector of means and variances.
Using the notation from Rasmussen & Williams, see e.g. equation 5.12:
    σᵢ = 𝕍 (yᵢ | y₋ᵢ)^(1/2)
    μᵢ = 𝔼 (yᵢ | y₋ᵢ)

The output is the same as fitting the GP on x₋ᵢ,y₋ᵢ,  and calling `predict_f`
on xᵢ to obtain the LOO predictive mean and variance, repeated for each observation i.
With GPs, this can thankfully be done analytically with a bit of linear algebra,
which is what this function implements.

See also: [`predict_CVfold`](@ref), [`logp_LOO`](@ref)
"""
function predict_LOO(gp::GPE)
    # extract relevant bits from GPE object
    Σ = gp.cK
    alpha = gp.alpha
    y = gp.y
    return predict_LOO(Σ, alpha, y)
end

"""
    logp_LOO(gp::GPE)

Leave-one-out log probability CV criterion.
This is implemented by  summing the
normal log-pdf of each observation
with predictive LOO mean and variance parameters
obtained by the `predict_LOO` function.

See also: [`logp_CVfold`](@ref), [`update_mll!`](@ref)
"""
function logp_LOO(gp::GPE)
    y = gp.y
    μ, σ2 = predict_LOO(gp)
    return sum(logpdf(Normal(μi,√σi2), yi)
               for (μi,σi2,yi)
               in zip(μ,σ2,y)
               )
end

"""
    dlogpdθ_LOO_kern!(gp::GPE)

Gradient of leave-one-out CV criterion with respect to the kernel hyperparameters.
See Rasmussen & Williams equations 5.13.

See also: [`logp_LOO`](@ref), [`dlogpdσ2_LOO`](@ref), [`dlogpdθ_LOO`](@ref)
"""
function dlogpdθ_LOO_kern!(∂logp∂θ::AbstractVector{<:Real}, invΣ::PDMat, kernel::Kernel, x::AbstractMatrix, y::AbstractVector, data::KernelData, alpha::AbstractVector)
    dim = num_params(kernel)
    nobs = length(y)
    @assert length(∂logp∂θ) == dim

    σi2 = 1 ./ diag(invΣ)
    μi = -alpha .* σi2 .+ y

    # Note: if useful, the derivatives of μ and σ could be moved to a separate function.
    # ∂μ∂θ = Matrix{Float64}(undef, nobs, dim)
    # ∂σ∂θ = Matrix{Float64}(undef, nobs, dim)
    Zj = Matrix{Float64}(undef, nobs, nobs)
    for j in 1:dim
        grad_slice!(Zj, kernel, x, x, data, j)
        Zj = invΣ.mat * Zj
        # ldiv!(Σ, Zj)

        ZjΣinv = diag(Zj*Matrix(invΣ))
        ∂σ2∂θj = ZjΣinv.*(σi2.^2)
        ∂μ∂θj = (Zj*alpha).*σi2 .- alpha .* ∂σ2∂θj
        # ∂μ∂θ[:,j] = ∂μ∂θj
        # ∂σ∂θ[:,j] = ∂σ∂θj

        ∂logp∂θj = 0.0
        for i in 1:nobs
            # exponentiated quadratic component:
            ∂logp∂θj -= 2*(y[i]-μi[i]) / σi2[i] * ∂μ∂θj[i]
            ∂logp∂θj -= (y[i]-μi[i])^2 * ZjΣinv[i]
            # log determinant component:
            @assert ZjΣinv[i] * σi2[i] ≈ ∂σ2∂θj[i] / σi2[i]
            ∂logp∂θj += ZjΣinv[i] * σi2[i]
        end
        ∂logp∂θ[j] = ∂logp∂θj
    end
    ∂logp∂θ .*= -1/2
    return ∂logp∂θ
end

"""
    dlogpdσ2_LOO(invΣ::PDMat, x::AbstractMatrix, y::AbstractVector, data::KernelData, alpha::AbstractVector)

Derivative of leave-one-out CV criterion with respect to the logNoise parameter.

See also: [`logp_LOO`](@ref), [`dlogpdθ_LOO_kern!`](@ref), [`dlogpdθ_LOO`](@ref)
"""
function dlogpdσ2_LOO(invΣ::PDMat, x::AbstractMatrix, y::AbstractVector, data::KernelData, alpha::AbstractVector)
    nobs = length(y)

    σi2 = 1 ./ diag(invΣ)
    μi = -alpha .* σi2 .+ y

    Zj = invΣ.mat
    ZjΣinv = diag(Zj^2)
    ∂σ2∂σ2 = ZjΣinv.*(σi2.^2)
    ∂μ∂σ2 = (Zj*alpha).*σi2 .- alpha .* ∂σ2∂σ2

    ∂logp∂σ2 = 0.0
    for i in 1:nobs
        # exponentiated quadratic component:
        ∂logp∂σ2 -= 2*(y[i]-μi[i]) / σi2[i] * ∂μ∂σ2[i]
        ∂logp∂σ2 -= (y[i]-μi[i])^2 * ZjΣinv[i]
        # log determinant component:
        @assert ZjΣinv[i] * σi2[i] ≈ ∂σ2∂σ2[i] / σi2[i]
        ∂logp∂σ2 += ZjΣinv[i] * σi2[i]
    end
    return -∂logp∂σ2 ./ 2
end

"""
    dlogpdθ_LOO(gp::GPE; noise::Bool, domean::Bool, kern::Bool)

Derivatives of the leave-one-out CV criterion.

See also: [`logp_LOO`](@ref), [`dlogpdθ_LOO_kern!`](@ref), [`dlogpdσ2_LOO`](@ref)
"""
function dlogpdθ_LOO(gp::GPE; noise::Bool, domean::Bool, kern::Bool)
    Σ = gp.cK
    x, y = gp.x, gp.y
    data = gp.data
    kernel = gp.kernel
    alpha = gp.alpha

    invΣ = inv(Σ)

    n_mean_params = num_params(gp.mean)
    n_kern_params = num_params(gp.kernel)
    ∂logp∂θ = Vector{Float64}(undef, noise + domean*n_mean_params + kern*n_kern_params)
    i = 1
    if noise
        ∂logp∂θ[i] = dlogpdσ2_LOO(invΣ, x, y, data, alpha)*2*exp(2 * gp.logNoise.value)
        i += 1
    end
    if domean && n_mean_params>0
        throw("I don't know how to do means yet")
        Mgrads = grad_stack(gp.mean, gp.x)
        for j in 1:n_mean_params
            gp.dmll[i] = dot(Mgrads[:,j], gp.alpha)
            i += 1
        end
    end
    if kern
        ∂logp∂θ_k = @view(∂logp∂θ[i:end])
        dlogpdθ_LOO_kern!(∂logp∂θ_k, invΣ, kernel, x, y, data, alpha)
    end
    return ∂logp∂θ
end

########################
# Arbitrary fold       #
########################

const Folds = AbstractVector{<:AbstractVector{Int}}

function predict_CVfold(Σ::AbstractPDMat, alpha::AbstractVector{<:Real}, y::AbstractVector{<:Real}, folds::Folds)
    invΣ = inv(Σ)
    μ = Vector{Float64}[]
    Σ = Matrix{Float64}[]
    for V in folds
        ΣVT = inv(Matrix(invΣ)[V,V])
        μVT = y[V]-ΣVT*alpha[V]
        push!(μ, μVT)
        push!(Σ, ΣVT)
    end
    return μ, Σ
end

"""
    predict_CVfold(gp::GPE, folds::Folds)

Cross-validated predictions for arbitrary folds.
A fold is a set of indices of the validation set.
Returns predictions of y_V (validation set) given all other observations y_T (training set),
as a vector of means and covariances.
Using the notation from Rasmussen & Williams, see e.g. equation 5.12:
    σᵢ = 𝕍 (yᵢ | y₋ᵢ)^(1/2)
    μᵢ = 𝔼 (yᵢ | y₋ᵢ)

The output is the same as fitting the GP on x_T,y_T,  and calling `predict_f`
on x_V to obtain the LOO predictive mean and covariance, repeated for each fold V.
With GPs, this can thankfully be done analytically with a bit of linear algebra,
which is what this function implements.

See also: [`predict_LOO`](@ref)
"""
function predict_CVfold(gp::GPE, folds::Folds)
    # extract relevant bits from GPE object
    Σ = gp.cK
    alpha = gp.alpha
    y = gp.y
    return predict_CVfold(Σ, alpha, y, folds)
end

"""
    logp_CVfold(gp::GPE, folds::Folds)

CV criterion for arbitrary fold.  A fold is a set of indices of the validation set.

See also [`predict_CVfold`](@ref), [`logp_LOO`](@ref)
"""
function logp_CVfold(gp::GPE, folds::Folds)
    y = gp.y
    μ, Σ = predict_CVfold(gp, folds)
    CV = 0.0
    for (μVT,ΣVT,V) in zip(μ,Σ,folds)
        chol = similar(ΣVT)
        ΣVT, chol = make_posdef!(ΣVT, chol; nugget=1e-10)
        ΣPD = PDMat(ΣVT, chol)
        CV += logpdf(MvNormal(μVT, ΣPD), y[V])
    end
    return CV
end

"""
    gradient_fold(invΣ, alpha, ZjΣinv, Zjα, V::AbstractVector{Int})

Gradient with respect to the kernel hyperparameters of the CV criterion
component for one validation fold:
```math
    \\nabla_{\\theta}\\left(\\log p (Y_V \\mid Y_T, \\theta) \\right)
```
where `Y_V` is the validation set and `Y_T` is the training set (all other observations)
for this fold.
"""
@inline function gradient_fold(invΣ, alpha, ZjΣinv, Zjα, V::AbstractVector{Int})
    ∂logp∂θj = 0.0
    ΣVTinv = PDMats.PDMat(mat(invΣ)[V,V])
    ΣVTα = ΣVTinv\alpha[V]
    # exponentiated quadratic component:
    ZjΣinvVV = ZjΣinv[V,V]
    ∂logp∂θj -= 2*dot(ΣVTα,Zjα[V])
    ∂logp∂θj += dot(ΣVTα, ZjΣinvVV*ΣVTα)
    # log determinant component:
    ∂logp∂θj += tr(ΣVTinv\ZjΣinvVV)
    return ∂logp∂θj
end

function dlogpdθ_CVfold_kern!(∂logp∂θ::AbstractVector{<:Real}, invΣ::PDMat, kernel::Kernel, x::AbstractMatrix, y::AbstractVector, data::KernelData, alpha::AbstractVector, folds::Folds)
    nobs = length(y)
    dim = num_params(kernel)

    @assert length(∂logp∂θ) == dim
    buffer1 = Matrix{Float64}(undef, nobs, nobs)
    buffer2 = Matrix{Float64}(undef, nobs, nobs)
    for j in 1:dim
        grad_slice!(buffer2, kernel, x, x, data, j)
        mul!(buffer1, mat(invΣ), buffer2)
        Zj = buffer1
        Zjα = Zj*alpha

        mul!(buffer2, Zj, mat(invΣ))
        ZjΣinv = buffer2

        ∂logp∂θj = 0.0
        for V in folds
            ∂logp∂θj += gradient_fold(invΣ, alpha, ZjΣinv, Zjα, V)
        end
        ∂logp∂θ[j] = ∂logp∂θj
    end
    ∂logp∂θ .*= -1/2
    return ∂logp∂θ
end

function dlogpdσ2_CVfold(invΣ::PDMat, x::AbstractMatrix, y::AbstractVector, data::KernelData, alpha::AbstractVector, folds::Folds)
    nobs = length(y)

    Zj = mat(invΣ)
    Zjα = Zj*alpha
    ZjΣinv = Zj^2

    ∂logp∂σ2 = 0.0
    for V in folds
        ∂logp∂σ2 += gradient_fold(invΣ, alpha, ZjΣinv, Zjα, V)
    end
    return -∂logp∂σ2 / 2
end

"""
    dlogpdθ_CVfold_kern!(∂logp∂θ::AbstractVector{<:Real}, gp::GPE, folds::Folds;
                         noise::Bool, domean::Bool, kern::Bool)

Derivative of leave-one-out CV criterion with respect to the noise, mean and kernel hyperparameters.
See Rasmussen & Williams equations 5.13.
"""
function dlogpdθ_CVfold(gp::GPE, folds::Folds; noise::Bool, domean::Bool, kern::Bool)
    Σ = gp.cK
    x, y = gp.x, gp.y
    data = gp.data
    kernel = gp.kernel
    alpha = gp.alpha

    invΣ = inv(Σ)

    n_mean_params = num_params(gp.mean)
    n_kern_params = num_params(gp.kernel)
    ∂logp∂θ = Vector{Float64}(undef, noise + domean*n_mean_params + kern*n_kern_params)
    i = 1
    if noise
        ∂logp∂θ[i] = dlogpdσ2_CVfold(invΣ, x, y, data, alpha, folds)*2*exp(2 * gp.logNoise.value)
        i += 1
    end
    if domean && n_mean_params>0
        throw("I don't know how to do means yet")
        Mgrads = grad_stack(gp.mean, gp.x)
        for j in 1:n_mean_params
            gp.dmll[i] = dot(Mgrads[:,j], gp.alpha)
            i += 1
        end
    end
    if kern
        ∂logp∂θ_k = @view(∂logp∂θ[i:end])
        dlogpdθ_CVfold_kern!(∂logp∂θ_k, invΣ, kernel, x, y, data, alpha, folds)
    end
    return ∂logp∂θ
end


text:
const Comprehension = Union{
    Iterators.Generator,
    Iterators.Filter,
    Iterators.Flatten,
}

"""
    eduction(iterator::Iterators.Generator)
    eduction(iterator::Iterators.Filter)
    eduction(iterator::Iterators.Flatten)

Convert an `iterator` to an eduction.  The iterators that are
typically used in the generator comprehensions are supported.

!!! compat "Transducers.jl 0.3"

    New in version 0.3.

# Examples
```jldoctest
julia> using Transducers

julia> iter = (y for x in 1:10 if x % 2 == 0 for y in (x, x + 1));

julia> ed = eduction(iter);

julia> collect(iter) == collect(ed)
true
```
"""
eduction(itr::Comprehension) = eduction(extract_transducer(itr)...)
eduction(ed::Eduction) = ed

"""
    Transducer(iterator::Iterators.Generator)
    Transducer(iterator::Iterators.Filter)
    Transducer(iterator::Iterators.Flatten)

Extract "processing" part of an `iterator` as a `Transducer`.  The
"data source" iterator (i.e., `xs` in `(f(x) for x in xs)`) is ignored
and `nothing` must be used as a place holder (i.e., `(f(x) for x in
nothing)`).

See also [`eduction`](@ref).

!!! compat "Transducers.jl 0.3"

    New in version 0.3.

# Examples
```jldoctest
julia> using Transducers

julia> xf1 = Transducer(2x for x in nothing if x % 2 == 0);

julia> xf2 = Filter(x -> x % 2 == 0) |> Map(x -> 2x);  # equivalent

julia> xs = 1:10
       collect(xf1, xs) == collect(xf2, xs)
true
```
"""
function Transducer(iter::Comprehension)
    xf, dataiter = extract_transducer(iter)
    dataiter === nothing && return xf

    if dataiter isa Iterators.ProductIterator
        throw(ArgumentError("""
              Conversion of the product generator comprehension, i.e.,
                  Transducer(... for x in xs, y in ys ...)
              is not supported.
              """))
    else
        throw(ArgumentError("""
              Unsupported iterator:
                  $(typeof(dataiter))
              Note that the "data source" iterator must be set to `nothing`.
              """))
    end
end

iterinner(iter::Iterators.Generator) = iter.iter
iterinner(iter::Iterators.Filter) = iter.itr
iterinner(iter::Iterators.Flatten) = iter.it
# iterf(iter::Iterators.Generate) = iter.f
# iterf(iter::Iterators.Filter) = iter.flt

"""
    extract_transducer(foldable) -> (xf, foldable′)

"Reverse" of [`eduction`](@ref).

# Examples
```jldoctest
julia> using Transducers

julia> double(x) = 2x;

julia> xs = 1:10;

julia> xf, foldable = Transducers.extract_transducer(double(x) for x in xs);

julia> xf == Map(double)
true

julia> foldable == xs
true
```
"""
extract_transducer(iter) = IdentityTransducer(), iter
extract_transducer(ed::Eduction) = Transducer(ed.rf), ed.coll

function extract_transducer(iter::Iterators.Generator)
    xf, bottom = extract_transducer(iterinner(iter))
    return xf |> Map(iter.f), bottom
end

function extract_transducer(iter::Iterators.Filter)
    xf, bottom = extract_transducer(iterinner(iter))
    return xf |> Filter(iter.flt), bottom
end

function extract_transducer(iter::Iterators.Flatten)
    xf, bottom = extract_transducer(iterinner(iter))
    return xf |> Cat(), bottom
end

"""
    Transducers.NoAdjoint(itr)

Bypass the optimization step by [`retransform`](@ref).
"""
struct NoAdjoint{T}
    itr::T
end

SplittablesBase.amount(itr::NoAdjoint) = amount(itr.itr)

extract_transducer(itr::NoAdjoint) = IdentityTransducer(), itr.itr

"""
    Transducers.retransform(rf, itr) -> rf′, itr′

Extract transformations in `rf` and `itr` and use the appropriate adjoint for
better performance.

# Examples
```jldoctest
julia> using Transducers

julia> double(x) = 2x;

julia> itr0 = 1:10;

julia> itr1 = (double(x) for x in itr0);

julia> rf, itr2 = Transducers.retransform(+, itr1);

julia> itr2 === itr0
true

julia> rf == reducingfunction(Map(double), +)
true
```
"""
function retransform(rf, itr1)
    xf, itr0 = extract_transducer(itr1)
    return reducingfunction(xf, rf), itr0
end

# TODO: Consider doing the "opposite" of `extract_transducer` for
# `Partition` etc.:
# https://github.com/JuliaFolds/Transducers.jl/issues/7


text:
using CSV
using DataFrames
using DelimitedFiles
using YAML
using Query

"""
    partition!(d::Dict, set::Dict; kwargs...)
"""
function partition!(d::Dict, set::Dict; save = true, overwrite = false)

    d_read = read_build("partition"; save = save, overwrite = overwrite);
    if !isempty(d_read)
        [d[k] = v for (k,v) in d_read]
        return d
    end

    [d[k] = edit_with(filter_with(d[k], (yr = set[:yr],)), [Drop(:units,"all","==")])
        for k in [:supply, :use]]

    _partition_io!(d, set)
    _partition_fd0!(d, set)
    _partition_ts0!(d, set)
    _partition_va0!(d, set)
    _partition_x0!(d, set)

    _partition_cif0!(d, set)
    _partition_duty0!(d, set)
    _partition_mrg0!(d, set)
    _partition_sbd0!(d, set)
    _partition_tax0!(d, set)

    _partition_m0!(d, set)   # cif0
    _partition_trn0!(d, set) # cif0

    _partition_fs0!(d)       # fd0
    _partition_s0!(d)        # ys0
    _partition_md0!(d, set)  # mrg0, trn0
    _partition_ms0!(d)       # mrg0, trn0

    _partition_y0!(d, set)   # ms0, fs0, ys0
    _partition_a0!(d, set)   # fd0, id0

    _partition_ta0!(d)       # a0, sbd0, tax0
    _partition_tm0!(d)       # duty0, m0

    # _partition_lshr0!(d, set) # va0

    d_save = delete!(delete!(copy(d), :supply), :use)
    write_build("partition", d_save; save = save)

    return d
end

"""
    _remove_imrg(df::DataFrame, x::Pair{Symbol,Array{String,1}})
Removes commodity taxes and subsidies on the goods which are produced solely
for supplying retail sales margin.
"""
function _remove_imrg(df::DataFrame, x::Pair{Symbol,Array{String,1}})
    df[findall(in(x.second), df[:,x.first]), :value] .= 0.0
    return df
end

"""
    _partition_io!(d::Dict, set::Dict)
"""
function _partition_io!(d::Dict, set::Dict)
    println("  Partitioning id0 and ys0, supply/demand data.")
    d[:id0] = filter_with(d[:use], set)
    d[:ys0] = filter_with(d[:supply], set)

    (d[:id0], d[:ys0]) = fill_zero(d[:id0], d[:ys0])

    # Treat negative inputs as outputs.
    d[:ys0][!,:value] = d[:ys0][:,:value] - min.(0, d[:id0][:,:value])
    d[:id0][!,:value] = max.(0, d[:id0][:,:value])

    d[:id0] = dropzero(d[:id0])
    d[:ys0] = dropzero(d[:ys0])
end

"""
    _partition_a0!(d::Dict, set::Dict)
`a0`: Armington supply
"""
function _partition_a0!(d::Dict, set::Dict)
    println("  Partitioning a0, Armington supply")
    d[:a0] = combine_over(d[:fd0], :fd) + combine_over(d[:id0], :j)
    d[:a0] = _remove_imrg(d[:a0], :i => set[:imrg])
end

"""
    _partition_bopdef0!(d::Dict, set::Dict)
`bopdef0`: Balance of payments deficit
"""
function _partition_bopdef0!(d::Dict, set::Dict)
    println("  Partitioning bopdef0, balance of payments deficit")
    d[:bopdef] = fill_zero((yr = set[:yr], ))
end

"""
    _partition_cif0!(d::Dict, set::Dict)
`cif0`: CIF/FOB Adjustments on Imports
"""
function _partition_cif0!(d::Dict, set::Dict)
    println("  Partitioning CIF/FOB adjustments on imports")
    d[:cif0] = filter_with(d[:supply], (i = set[:i], j = "ciffob"))[:,[:yr,:i,:value]]
end

"""
    _partition_duty0!(d::Dict, set::Dict)
`duty0`: Import duties
"""
function _partition_duty0!(d::Dict, set::Dict)
    println("  Partitioning duty0, import duties")
    d[:duty0] = filter_with(d[:supply], (i = set[:i], j = "duties"))[:,[:yr,:i,:value]]
    d[:duty0] = _remove_imrg(d[:duty0], :i => set[:imrg])
end

"""
    _partition_fd0!(d::Dict, set::Dict)
`fd0`: Final demand
"""
function _partition_fd0!(d::Dict, set::Dict)
    println("  Partitioning fd0, final demand")
    d[:fd0] = filter_with(d[:use], (i = set[:i],  j = set[:fd]))
    d[:fd0] = edit_with(d[:fd0], Rename(:j, :fd))
end

"""
    _partition_fs0!(d::Dict)
`fs0`: Household supply.
Move household supply of recycled goods into the domestic output market,
from which some may be exported.
"""
function _partition_fs0!(d::Dict)
    println("  Partitioning fs0, household supply")
    d[:fs0] = filter_with(d[:fd0], (fd = "pce",))[:,[:yr,:i,:value]]
    d[:fs0][!,:value] .= - min.(d[:fs0][:,:value], 0)
end

# """
#     _partition_lshr0!(d::Dict)
# `lshr0`: Labor share of value added
# """
# function _partition_lshr0!(d::Dict, set::Dict)
#     va0 = edit_with(unstack(copy(d[:va0]), :va, :value),
#         [Rename(:j,:s); Replace.(Symbol.(set[:va]), missing, 0.0); Drop(:units,"all","==")])

#     d[:lshr0]  = va0[:,[:yr,:s,:compen]]
#     d[:lshr0] /= (va0[:,[:yr,:s,:compen]] + va0[:,[:yr,:s,:surplus]])

#     # !!!!! _partition_lshr0 needs to come after calibration.
#     # Order is: io, calibrate, share, disagg.
#     d[:lshr0][va0[:,:surplus] .< 0,:value] .== 1.0
#     dropmissing!(d[:lshr0])
# end

"""
    _partition_m0!(d::Dict, set::Dict)
`m0`: Imports
"""
function _partition_m0!(d::Dict, set::Dict)
    println("  Partitioning m0, imports")
    d[:m0]  = filter_with(d[:supply], (i = set[:i], j = "imports"))[:,[:yr,:i,:value]]

    # Adjust transport margins for transport sectors according to CIF/FOB adjustments.
    # Insurance imports are specified as net of adjustments.
    d[:m0] += filter_with(d[:cif0], (i = "ins",))
    d[:m0] = _remove_imrg(d[:m0], :i => set[:imrg])
end

"""
    _partition_md0!(d::Dict)
`md0`: Margin demand
"""
function _partition_md0!(d::Dict, set::Dict)
    println("  Partitioning md0, margin demand")
    d[:md0] = [edit_with(d[:mrg0], Add(:m, "trd")); edit_with(d[:trn0], Add(:m, "trn"))]
    d[:md0] = sort(d[:md0][:,[:yr,:i,:m,:value]])
    d[:md0] = _remove_imrg(d[:md0], :i => set[:imrg])

    d[:md0][!,:value] .= max.(d[:md0][:,:value], 0)
end

"""
    _partition_ms0!(d::Dict)
`ms0`: Margin supply
"""
function _partition_ms0!(d::Dict)
    println("  Partitioning ms0, margin supply")
    d[:ms0] = [edit_with(d[:mrg0], Add(:m, "trd")); edit_with(d[:trn0], Add(:m, "trn"))]
    d[:ms0] = sort(d[:ms0][:,[:yr,:i,:m,:value]])

    d[:ms0][!,:value] .= max.(-d[:ms0][:,:value], 0)
end

"""
    _partition_mrg0!(d::Dict, set::Dict)
`mrg0`: Trade margins
"""
function _partition_mrg0!(d::Dict, set::Dict)
    println("  Partitioning mrg0, trade margins")
    d[:mrg0] = filter_with(d[:supply], (i = set[:i], j = "margins"))[:,[:yr,:i,:value]]
end

"""
    _partition_s0!(d::Dict)
`s0`: Aggregate supply
"""
function _partition_s0!(d::Dict)
    println("  Partitioning s0, aggregate supply")
    d[:s0] = combine_over(d[:ys0], :i)
end

"""
    _partition_sbd0!(d::Dict, set::Dict)
`sbd0`: Subsidies on products
"""
function _partition_sbd0!(d::Dict, set::Dict)
    println("  Partitioning sbd0, subsidies on products")
    d[:sbd0] = filter_with(d[:supply], (i = set[:i], j = "subsidies"))[:,[:yr,:i,:value]]
    d[:sbd0] = _remove_imrg(d[:sbd0], :i => set[:imrg])
    d[:sbd0][!,:value] *= -1
end

"""
    _partition_ta0!(d::Dict)
`ta0`: Import tariff
"""
function _partition_ta0!(d::Dict)
    println("  Partitioning ta0, import tariffs")
    d[:ta0] = dropnan((d[:tax0] - d[:sbd0]) / d[:a0])
    # d[:ta0] = edit_with(d[:ta0], Drop(:units,"all","=="))
end

"""
    _partition_tax0!(d::Dict, set::Dict)
`tax0`: Taxes on products
"""
function _partition_tax0!(d::Dict, set::Dict)
    println("  Partitioning tax0, taxes on products")
    d[:tax0] = filter_with(d[:supply], (i = set[:i], j = "tax"))[:,[:yr,:i,:value]]
    d[:tax0] = _remove_imrg(d[:tax0], :i => set[:imrg])
end

"""
    _partition_tm0!(d::Dict)
`tm0`: Tax net subsidy rate on intermediate demand
"""
function _partition_tm0!(d::Dict)
    println("  Partitioning tm0, tax net subsidy rate on intermediate demand")
    d[:tm0] = dropnan(d[:duty0] / d[:m0])
    # d[:tm0] = edit_with(d[:tm0], Drop(:units,"all","=="))
end

"""
    _partition_trn0!(d::Dict, set::Dict)
`trn0`: Transportation costs
"""
function _partition_trn0!(d::Dict, set::Dict)
    println("  Partitioning trn0, transportation costs")
    d[:trn0]  = filter_with(d[:supply], (i = set[:i], j = "trncost"))[:,[:yr,:i,:value]]

    # Adjust transport margins for transport sectors according to CIF/FOB adjustments.
    # Insurance imports are specified as net of adjustments.
    d[:trn0] += edit_with(d[:cif0], Drop(:i,"ins","=="))
end

"""
    _partition_ts0!(d::Dict, set::Dict)
`ts0`: Taxes and subsidies
"""
function _partition_ts0!(d::Dict, set::Dict)
    println("  Partitioning ts0, taxes and subsidies")
    d[:ts0] = filter_with(d[:use], (i = set[:ts], j = set[:j]))
    d[:ts0][d[:ts0][:,:i] .== "subsidies", :value] *= -1  # treat negative inputs as outputs    return d
end

"""
    _partition_va0!(d::Dict, set::Dict)
`va0`: Value added
"""
function _partition_va0!(d::Dict, set::Dict)
    println("  Partitioning va0, value added")
    d[:va0] = filter_with(d[:use], (i = set[:va], j = set[:j]))
    d[:va0] = edit_with(d[:va0], Rename(:i, :va))
end

"""
    _partition_x0!(d::Dict, set::Dict)
`x0`: Exports of goods and services
"""
function _partition_x0!(d::Dict, set::Dict)
    println("  Partitioning x0, exports of goods and services")
    d[:x0] = filter_with(d[:use], (i = set[:i], j = "exports"))[:,[:yr,:i,:value]]
    d[:x0] = _remove_imrg(d[:x0], :i => set[:imrg])
end

"""
    _partition_y0!(d::Dict, set::Dict)
`y0`: Gross output
"""
function _partition_y0!(d::Dict, set::Dict)
    println("  Partitioning y0, gross output")
    d[:y0] = combine_over(d[:ys0], :j) + d[:fs0] - combine_over(d[:ms0], :m)
    d[:y0] = _remove_imrg(d[:y0], :i => set[:imrg])
end

max_stars_repo_name: matbesancon/SLiDE

max_issues_repo_name: matbesancon/SLiDE

max_forks_repo_name: matbesancon/SLiDE


text:
using Statistics

function show_full(array)
    show(IOContext(stdout, :limit=>false), MIME"text/plain"(), array)
end

## basic simulation functions
function xdot(x,y,z; sigma=10)
    return sigma*(y-x)
end
function ydot(x,y,z; rho=28)
    return x*(rho-z)-y
end
function zdot(x,y,z; beta=8/3)
    return x*y - beta*z
end

function make_step(x,y,z; delta=0.01)
    xp = x + xdot(x,y,z)*delta
    yp = y + ydot(x,y,z)*delta
    zp = z + zdot(x,y,z)*delta

    xn = x + 0.5*(xdot(x,y,z)+xdot(xp,yp,zp))*delta
    yn = y + 0.5*(ydot(x,y,z)+ydot(xp,yp,zp))*delta
    zn = z + 0.5*(zdot(x,y,z)+zdot(xp,yp,zp))*delta

    return xn, yn, zn
end

## threshold functions

function threshold_basic(x,y,z)
    # Z=24
    out = 0
    if z>24
        out = 1
    end
    return out
end

function z_tercile_thresh(x,y,z)
    # states should be mutually exclusive
    # Z=18.3 and Z=32
    if z<18.3 #bottom
        return 1
    end
    if z<32 #middle
        return 2
    end
    return 3 #top
end

function eight_state_thresh(x,y,z)
    state = [0,0,0]
    if x > 0
        state[1] = 1
    end
    if y > 0
        state[2] = 1
    end
    if z > 24
        state[3] = 1
    end
    map = Dict([0,0,0] => 1, [1,0,0] => 2, [0,1,0] => 3, [1,1,0] => 4,
                [0,0,1] => 5, [1,0,1] => 6, [0,1,1] => 7, [1,1,1] => 8)
    return  map[state]
end

function twelve_state_quant_thresh(x,y,z)
    #using 0.05, 0.95 quantiles for z
    state = [0,0,0]
    if x > 0
        state[1] = 1
    end
    if y > 0
        state[2] = 1
    end
    if z > 37.9
        state[3] = 2
    elseif z > 10.6
        state[3] = 1
    end
    # inefficient to be defining the map inside the function...
    map = Dict([0,0,0] => 1, [1,0,0] => 2, [0,1,0] => 3, [1,1,0] => 4,
    [0,0,1] => 5, [1,0,1] => 6, [0,1,1] => 7, [1,1,1] => 8,
    [0,0,2] => 9, [1,0,2] => 10, [0,1,2] => 11, [1,1,2] => 12)
    return map[state]
end

## simulation run
function run_sim(;runs=3000, timing=false, thresh_func=z_tercile_thresh)
    x, y, z = 0, 1, 0
    X, Y, Z = Float64[0,], Float64[1,], Float64[0,]
    cnt = 0
    state_list = [] #list of unique states
    holding_times = [] #list of holding times of each unique state

    state = thresh_func(x,y,z)

    for i in 1:runs
        x,y,z = make_step(x,y,z)
        push!(X, x)
        push!(Y, y)
        push!(Z, z)

        if thresh_func(x,y,z) == state
            cnt += 1
        else
            push!(state_list, state)
            push!(holding_times, cnt)
            state = thresh_func(x,y,z)
            cnt = 0
        end

    end
    return X, Y, Z, state_list, holding_times
end

X, Y, Z, state_list, holding_times = run_sim(;runs=1000000, timing=true, thresh_func=twelve_state_quant_thresh)
number_of_states = 12

#basic plotting
using Plots
Plots.plot(X,Y)
Plots.plot(Z[100:end])
Plots.histogram(Z;bins=100)

#construct probability transition matrix directly
function construct_p(number_of_states, state_list, holding_times)
    constr_prob = zeros(number_of_states, number_of_states)

    for i in 1:length(state_list)-1
        local current = state_list[i]
        local next = state_list[i+1]
        constr_prob[current, current] += holding_times[i]-1
        constr_prob[next, current] += 1 #reconstruct full markov chain
    end
    norm = sum(constr_prob, dims=1)

    norm_constr_prob = constr_prob ./ norm
    empirical_T = log(norm_constr_prob)/0.01
    return (norm_constr_prob, empirical_T)
end
norm_constr_prob, empirical_T = construct_p(number_of_states, state_list, holding_times)
show_full(norm_constr_prob)
show_full(empirical_T)

## other method!
function construct_t(number_of_states, state_list, holding_times)
    holding_time_dist = [[] for n in 1:number_of_states]
    constr_T = zeros(number_of_states, number_of_states)

    #count number of times unique states appear
    for i in 1:length(state_list)-1 #for i in unique states #why is the -1??
        local current = state_list[i]
        local next = state_list[i+1]
        constr_T[next, current] += 1
        push!(holding_time_dist[current], holding_times[i]*0.01)
    end

    norm = sum(constr_T, dims=1)
    norm_constr_T = constr_T ./ norm

    holding_scale = 1 ./ mean.(holding_time_dist)
    for i in 1:number_of_states
        norm_constr_T[i, i] = -1.0 #this is the bit I don't quite understand
        norm_constr_T[:, i] *= holding_scale[i]
    end

    empirical_prob = exp(norm_constr_T * 0.01)
    return (norm_constr_T, empirical_prob)
end

norm_constr_T, empirical_prob = construct_t(number_of_states, state_list, holding_times)
show_full(norm_constr_T)
show_full(empirical_prob)

`                                                           `

## identify quantiles in Z
z_sort = copy(Z)
three_state_extremes = quantile!(z_sort,[0.05, 0.95]) #output: [10.6, 37.9]
five_state_extremes = quantile!(z_sort,[0.05, 0.25, 0.75, 0.95]) #output: [10.6, 16.9, 30.6, 37.9]

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

function get_maxima(ls; add_min=false)
    maxs = Float64[]
    if add_min
        mins = Float64[]
        if ls[1]<ls[2]
            push!(mins,ls[1])
        end
    end
    if ls[1]>ls[2]
        push!(maxs,ls[1])
    end
    for i=2:length(ls)-1
        if ls[i-1]<ls[i]>ls[i+1]
            push!(maxs,ls[i])
        end
        if add_min
            if ls[i-1]>ls[i]<ls[i+1]
                push!(mins,ls[i])
            end
        end
    end
    if ls[end]>ls[end-1]
        push!(maxs,ls[end])
    end
    if add_min
        if ls[end]<ls[end-1]
            push!(mins,ls[end])
        end
        return maxs, mins
    else
        return maxs
    end
end

# get subsequent maxima  "tent" plot
z_max = get_maxima(Z)
z_plot = Tuple{Float64, Float64}[]
for i=1:length(z_max)-1
    push!(z_plot,(z_max[i], z_max[i+1]))
end
i1 = [x[1] for x in z_plot]
i2 = [x[2] for x in z_plot]
Plots.scatter(i1,i2)

# get sequential max/min plot
z_max, z_min = get_maxima(Z, add_min=true)
Plots.scatter(z_max, z_min)

max_stars_repo_name: m-geo/lorenz

max_issues_repo_name: m-geo/lorenz

max_forks_repo_name: m-geo/lorenz


text:
# the following provides a deferred access mechanism to handle nested aggregate fields (in aggregates or arrays) to support correct/efficient behavior of:
#   a.b[3].c.d = x
#   y = a.b[3].c.d
const Cdeferrable = Union{Caggregate, Carray}
struct Caccessor{FieldType<:Union{Cdeferrable, Cconst{<:Cdeferrable}}, BaseType<:Union{Caggregate, Carray, Ptr{<:Caggregate}, Ptr{<:Carray}, Cconst{<:Caggregate}, Cconst{<:Carray}}, Offset<:Val}
	base::BaseType

	Caccessor{FieldType}(b::BaseType, ::Val{Offset} = Val(0)) where {FieldType, BaseType, Offset} = new{FieldType, BaseType, Val{Offset}}(b)
end

Base.convert(::Type{T}, ca::Caccessor{T}) where {T} = ca[]
Base.show(io::IO, ca::Caccessor) = show(io, ca[])

Base.getindex(cc::Caccessor{CC}) where {CC<:Cconst{<:Cdeferrable}} = Cconst{nonconst(CC)}(_bytes(cc))
Base.getindex(ca::Caccessor{CD}) where {CD<:Cdeferrable} = unsafe_load(_pointer(ca))
Base.setindex!(ca::Caccessor{CD}, val::CD) where {CD<:Cdeferrable} = unsafe_store!(_pointer(ca), val)

# Caggregate interface
const Caggregates = Union{CA, Cconst{CA}, Caccessor{CA}, Caccessor{<:Cconst{CA}}} where {CA<:Caggregate}
Base.propertynames(ca::CA; kwargs...) where {CA<:Caggregates} = propertynames(typeof(ca); kwargs...)
Base.propertynames(::Type{CA}; kwargs...) where {CA<:Caggregates} = map(((sym, fld),) -> sym, (sort(collect(Ctypelayout(_fieldtype(CA)).fields), by = ((sym, fld),) -> fld.ind)...,))

propertytypes(ca::CA; kwargs...) where {CA<:Caggregates} = propertytypes(typeof(ca); kwargs...)
propertytypes(::Type{CA}; kwargs...) where {CA<:Caggregates} = map(((sym, fld),) -> fld.type, (sort(collect(Ctypelayout(_fieldtype(CA)).fields), by = ((sym, fld),) -> fld.ind)...,))

Base.fieldnames(ca::CA; kwargs...) where {CA<:Caggregates} = fieldnames(typeof(ca); kwargs...)
Base.fieldnames(::Type{CA}; kwargs...) where {CA<:Caggregates} = propertynames(_fieldtype(CA); kwargs...)

Base.getproperty(ca::CA, sym::Symbol) where {CA<:Caggregates} = _getproperty(CA, _base(ca), Val(_fieldoffset(ca)), Ctypespec(_fieldtype(ca)), Val(sym))
Base.setproperty!(ca::CA, sym::Symbol, val) where {CA<:Caggregates} = _setproperty!(CA, _base(ca), Val(_fieldoffset(ca)), Ctypespec(_fieldtype(ca)), Val(sym), val)

# Carray interface
const Carrays = Union{CA, Cconst{CA}, Caccessor{CA}, Caccessor{<:Cconst{CA}}} where {CA<:Carray}
Base.getindex(ca::CA, ind) where {T<:Union{Cdeferrable, Cconst{<:Cdeferrable}}, N, _CA<:Carray{T, N}, CA<:Carrays{_CA}} = Caccessor{T}(_base(ca), Val(_fieldoffset(ca) + (ind-1)*sizeof(T)))
Base.getindex(ca::CA, ind) where {T, N, _CA<:Carray{T, N}, CA<:Carrays{_CA}} = unsafe_load(reinterpret(Ptr{T}, _pointer(ca)), ind)
Base.setindex!(ca::CA, val, ind) where {T, N, _CA<:Carray{T, N}, CA<:Carrays{_CA}} = unsafe_store!(reinterpret(Ptr{T}, _pointer(ca)), val, ind)

Base.firstindex(ca::CA) where {CA<:Carrays} = 1
Base.lastindex(ca::CA) where {CA<:Carrays} = length(ca)

Base.IndexStyle(::Type{CA}) where {CA<:Carrays} = IndexLinear()
Base.size(ca::CA) where {CA<:Carrays} = size(typeof(ca))
Base.length(ca::CA) where {CA<:Carrays} = length(typeof(ca))
Base.eltype(ca::CA) where {CA<:Carrays} = eltype(typeof(ca))
Base.size(::Type{CA}) where {T, N, _CA<:Carray{T, N}, CA<:Carrays{_CA}} = (N,)
Base.length(::Type{CA}) where {T, N, _CA<:Carray{T, N}, CA<:Carrays{_CA}} = N
Base.eltype(::Type{CA}) where {T, N, _CA<:Carray{T, N}, CA<:Carrays{_CA}} = T

Base.keys(ca::CA) where {CA<:Carrays} = firstindex(ca):lastindex(ca)
Base.values(ca::CA) where {CA<:Carrays} = iterate(ca)
Base.iterate(ca::CA, state = 1) where {CA<:Carrays} = state > length(ca) ? nothing : (ca[state], state+1)

_fieldoffset(cx::Union{Cdeferrable, Cconst, Caccessor}) = _fieldoffset(typeof(cx))
_fieldoffset(::Type{CD}) where {CD<:Cdeferrable} = 0
_fieldoffset(::Type{CC}) where {T, CC<:Cconst{T}} = _fieldoffset(T)
_fieldoffset(::Type{Caccessor{FieldType, BaseType, Val{Offset}}}) where {FieldType, BaseType, Offset} = Offset

_fieldtype(cx::Union{Cdeferrable, Cconst, Caccessor}) = _fieldtype(typeof(cx))
_fieldtype(::Type{CD}) where {CD<:Cdeferrable} = CD
_fieldtype(::Type{CC}) where {T, CC<:Cconst{T}} = _fieldtype(T)
_fieldtype(::Type{Caccessor{FieldType, BaseType, Val{Offset}}}) where {FieldType, BaseType, Offset} = FieldType

_base(cx::Union{Cdeferrable, Cconst}) = cx
_base(ca::Caccessor) = getfield(ca, :base)

_pointer(ptr::Ptr) = ptr
_pointer(cx::Union{Caggregate, Carray}) = reinterpret(Ptr{typeof(cx)}, pointer_from_objref(cx))
_pointer(ca::Caccessor) = reinterpret(Ptr{_fieldtype(ca)}, _pointer(_base(ca)) + _fieldoffset(ca))

_uint(::Type{T}) where {T} = sizeof(T) == sizeof(UInt8) ? UInt8 : sizeof(T) == sizeof(UInt16) ? UInt16 : sizeof(T) == sizeof(UInt32) ? UInt32 : sizeof(T) == sizeof(UInt64) ? UInt64 : sizeof(T) == sizeof(UInt128) ? UInt128 : error("Unable to create a UInt of $(sizeof(T)*8) bits")

function _bitmask(::Type{uint}, bits::Int) where {uint}
	mask = zero(uint)
	for i in 1:bits
		mask = (mask << one(uint)) | one(uint)
	end
	return uint(mask)
end

_readbyte(::Type{T}, base, ind) where {T<:Cconst} = :(getfield($(base), :mem)[$(ind)])
_readbyte(::Type{T}, base, ind) where {T<:Ptr} = :(unsafe_load(reinterpret(Ptr{UInt8}, $(base)), $(ind)))
_writebyte(::Type{T}, base, ind, val) where {T<:Ptr} = :(unsafe_store!(reinterpret(Ptr{UInt8}, $(base)), $(val), $(ind)))

_bytes(ca::Caccessor) = _bytes(_base(ca), Val(_fieldoffset(ca)), Val(sizeof(_fieldtype(ca))))
_bytes(cd::Cdeferrable, ::Val{offset}, ::Val{size}) where {offset, size} = _bytes(pointer_from_objref(cd), Val(offset), Val(size))
@generated function _bytes(base::Union{Cconst, Ptr}, ::Val{offset}, ::Val{size}) where {offset, size}
	return :($(map(ind -> _readbyte(base, :base, ind), offset+1:offset+size)...),)
end

@generated function _unsafe_load(base::Union{Cconst, Ptr}, ::Val{offset}, ::Type{uint}, ::Val{offbits}, ::Val{numbits}) where {offset, uint, offbits, numbits}
	sym = gensym("bitfield")
	result = [:($(sym) = uint(0))]
	for i in 1:sizeof(uint)
		todo"verify correctness on big endian machine"  #$((ENDIAN_BOM != 0x04030201 ? (sizeof(uint)-i) : (i-1))*8)
		offbits <= i*8 && (i-1)*8 < offbits+numbits && push!(result, :($(sym) |= uint($(_readbyte(base, :base, offset+i))) << uint($((i-1)*8))))
	end
	return quote let ; $(result...) ; $(sym) end end
end

@generated function _unsafe_store!(base::Ptr, ::Val{offset}, ::Type{uint}, ::Val{offbits}, ::Val{numbits}, val::uint) where {offset, uint, offbits, numbits}
	result = []
	for i in 1:sizeof(uint)
		todo"verify correctness on big endian machine"  #$((ENDIAN_BOM != 0x04030201 ? (sizeof(uint)-i) : (i-1))*8)
		offbits <= i*8 && (i-1)*8 < offbits+numbits && push!(result, _writebyte(base, :base, offset+i, :(UInt8((val >> $((i-1)*8)) & 0xff))))
	end
	return quote $(result...) end
end

@generated function _getproperty(::Type{<:Any}, base::Union{CA, Cconst{CA}, Ptr{CA}}, ::Val{offset}, ::Type{spec}, ::Val{name}) where {CA<:Caggregate, offset, spec<:Ctypespec, name}
	mem = base <: CA ? :(pointer_from_objref(base)) : :(base)

	fields = Ctypelayout(spec).fields
	if haskey(fields, name)
		typ = fields[name].type
		bits = fields[name].size
		off = fields[name].offset + offset*8

		if typ <: Union{Cdeferrable, Cconst{<:Cdeferrable}}
			return :(Caccessor{$(base <: Cconst || type(spec) <: Cconst ? Cconst(typ) : typ)}(base, Val($(off÷8))))
		else
			t = nonconst(typ)
			o = off & (8-1)
			b = iszero(bits) ? sizeof(t)*8 : bits
			uint = _uint(t)
			mask = _bitmask(uint, b)
			return quote
				val = _unsafe_load($(mem), Val($(off÷8)), $(uint), Val($(o)), Val($(b)))
				if $(!iszero(bits))
					val = (val >> $(uint)($(o))) & $(uint)($(mask))
					if $(t) <: Signed && ((val >> $(b-1)) & 1) != 0  # 0 = pos, 1 = neg
						val |= ~$(uint)(0) & ~$(uint)($(mask))
					end
				end
				return reinterpret($(t), val)
			end
		end
	end
	return :(error("Unable to get property `$(name)`, it is not a field of $(typeof(base))"))
end

_initproperty!(cx::CX, sym::Symbol, val) where {CA<:Caggregate, CX<:Union{CA, Caccessor{CA}}} = _setproperty!(CA, _base(cx), Val(_fieldoffset(cx)), Ctypespec(_fieldtype(cx)), Val(sym), val, Val(true))
@generated function _setproperty!(::Type{T}, base::Union{CA, Cconst{CA}, Ptr{CA}}, ::Val{offset}, ::Type{spec}, ::Val{name}, val, ::Val{force} = Val(false)) where {X, T<:Union{X, Caccessor{X}}, CA<:Caggregate, offset, spec<:Ctypespec, name, force}
	(X <: Cconst || base <: Cconst) && return :(error("Unable to change the value of a Cconst field"))
	mem = base <: CA ? :(pointer_from_objref(base)) : :(base)

	fields = Ctypelayout(spec).fields
	if haskey(fields, name)
		typ = fields[name].type
		bits = fields[name].size
		off = fields[name].offset + offset*8

		typ = force ? nonconst(typ) : typ
		typ <: Cconst && return :(error("Unable to change the value of a Cconst field"))

		if typ <: Carray
			return quote
				arr = Caccessor{$(typ)}(base, Val($(off÷8)))
				length(val) == length(arr) || error("Length of value does not match the length of the array field it is being assigned to")
				for (i, v) in enumerate(val)
					arr[i] = v
				end
				return val
			end
		elseif !iszero(bits)
			o = off & (8-1)
			b = iszero(bits) ? sizeof(typ)*8 : bits
			uint = _uint(typ)
			mask = _bitmask(uint, b) << o
			return quote
				field = _unsafe_load($(mem), Val($(off÷8)), $(uint), Val($(o)), Val($(b)))
				field &= ~$(uint)($(mask))
				field |= (reinterpret($(uint), convert($(typ), val)) << $(o)) & $(uint)($(mask))
				_unsafe_store!($(mem), Val($(off÷8)), $(uint), Val($(o)), Val($(b)), field)
				return val
			end
		else
			return quote
				unsafe_store!(reinterpret(Ptr{$(typ)}, $(mem) + $(off÷8)), val)
				return val
			end
		end
	end
	return :(error("Unable to set property `$(name)`, it is not a field of $(typeof(base))"))
end


text:
# Input data
data = map(x -> parse(Int64, x), split(readline("Day06.input"), ","));

# Part 1
getans =
    (days, data) -> sum(
        reduce(
            (acc, _) -> (v = circshift(acc, -1); v[7] += v[9]; v),
            1:days,
            init = map(x -> count(==(x - 1), data), 1:9),
        ),
    );

ans = getans(80, data);
println("Part 1: $ans")

# Part 2
ans = getans(256, data);
println("Part 2: $ans")


text:
using DelimitedFiles

inputs = Array{Int}(undef, 1);
outputs = [];

function modalize(code, parameters, modes)
    # an effort to make the ternary instructions a little less repetitive.
    left = modes[1] == 1 ? parameters[1] : code[parameters[1] + 1];
    right = modes[2] == 1 ? parameters[2] : code[parameters[2] + 1];
    return (left, right)
end

function intcodeAdd(code, parameters, modes, p::Int)
    (left, right) = modalize(code, parameters, modes)
    code[parameters[3] + 1] = left + right
    return p + intcode[code[p] % 100].n
end

function intcodeMultiply(code, parameters, modes, p::Int)
    (left, right) = modalize(code, parameters, modes)
    code[parameters[3] + 1] = left * right
    return p + intcode[code[p] % 100].n
end

function intcodeInput(code, parameters, modes, p::Int)
    code[parameters[1] + 1] = popfirst!(inputs)
    return p + intcode[code[p] % 100].n
end

function intcodeOutput(code, parameters, modes, p::Int)
    left = modes[1] == 1 ? parameters[1] : code[parameters[1] + 1];
    #println("Intcode output: $(left)")
    push!(outputs, left)
    return p + intcode[code[p] % 100].n
end

function intcodeJumpIfTrue(code, parameters, modes, p::Int)
    (left, right) = modalize(code, parameters, modes)
    if left != 0
        return right + 1
    else
        return p + intcode[code[p] % 100].n
    end
end

function intcodeJumpIfFalse(code, parameters, modes, p::Int)
    (left, right) = modalize(code, parameters, modes)
    if left == 0
        return right + 1
    else
        return p + intcode[code[p] % 100].n
    end
end

function intcodeLessThan(code, parameters, modes, p::Int)
    (left, right) = modalize(code, parameters, modes)
    code[parameters[3] + 1] = Int(left < right);
    return p + intcode[code[p] % 100].n
end

function intcodeEquals(code, parameters, modes, p::Int)
    (left, right) = modalize(code, parameters, modes)
    code[parameters[3] + 1] = Int(left == right);
    return p + intcode[code[p] % 100].n
end

function intcodeExit(code, parameters, modes, p::Int)
    Base.identity("Do nothing")
    return p + intcode[code[p] % 100].n
end

const intcode = Dict([
    (1, (f=intcodeAdd, n=4, name="Add")),
    (2, (f=intcodeMultiply, n=4, name="Multiply")),
    (3, (f=intcodeInput, n=2, name="Input")),
    (4, (f=intcodeOutput, n=2, name="Output")),
    (5, (f=intcodeJumpIfTrue, n=3, name="Jump-if-true")),
    (6, (f=intcodeJumpIfFalse, n=3, name="Jump-if-false")),
    (7, (f=intcodeLessThan, n=4, name="Less than")),
    (8, (f=intcodeEquals, n=4, name="Equals")),
    (99, (f=intcodeExit, n=1, name="Exit"))
]);

function run(code, noun=missing, verb=missing)
    c = copy(code); # copy so we can mutate the input
    c[2] = noun;
    c[3] = verb;
    return run(c);
end

function run(code)
    p::Int = 1; # instruction pointer
    global outputs = []
    while (inst = code[p]) != 99
        #println(c)
        opcode = inst % 100;
        parameters = view(code, (p + 1):(p + intcode[opcode].n - 1));
        modes = ((inst ÷ 100) % 10,
            (inst ÷ 1000) % 10,
            (inst ÷ 10000) % 10);
        #println("{IP=$(p)} $(intcode[opcode].name) ($(inst)): $(parameters)");
        p = intcode[opcode].f(code, parameters, modes, p);
    end
    #println(c)
    return code
end

const day5_test = ["1002,4,3,4,33"];
const day5_test_result = run(parse.(Int,split(day5_test[1],",")));
println("Day 5 Test: $(day5_test_result) ($(day5_test_result == [1002, 4, 3, 4, 99] ? "passed" : "failed"))");

inputs = [1];
input = readdlm(ARGS[1], ',', Int, '\n');
input = vec(input)
const day5_part1 = run(input);
println("Day 5 Part 1: $(pop!(outputs))");

const day5_part2_tests = ["3,9,8,9,10,9,4,9,99,-1,8", "3,9,7,9,10,9,4,9,99,-1,8",
    "3,3,1108,-1,8,3,4,3,99", "3,3,1107,-1,8,3,4,3,99"]
for test in day5_part2_tests
    global inputs = [1]
    global outputs = []
    result = run(parse.(Int,split(test,",")));
    println("Day 5 Test: $(result) produces outputs $(outputs)")
end

const day5_part2_jump = ["3,12,6,12,15,1,13,14,13,4,13,99,-1,0,1,9", "3,3,1105,-1,9,1101,0,0,12,4,12,99,1"];
for test in day5_part2_jump
    for i in 0:1
        global inputs = [i]
        result = run(parse.(Int,split(test,",")));
        println("Day 5 Jump Test: $(result) gave outputs $(outputs) for input=[$(i)]")
    end
end

const day5_large = "3,21,1008,21,8,20,1005,20,22,107,8,21,20,1006,20,31,
1106,0,36,98,0,0,1002,21,125,20,4,20,1105,1,46,104,
999,1105,1,46,1101,1000,1,20,4,20,1105,1,46,98,99"

for i in 7:9
    global inputs = [i];
    run(parse.(Int,split(day5_large, ",")))
    println("Day 5 Long Test: $(i): $(outputs)")
end

inputs = [5];
input = readdlm(ARGS[1], ',', Int, '\n');
input = vec(input)
run(input)
println("Day 5 Part 2: $(pop!(outputs))")


text:
int_rules_8_2 = @theory begin

    #= ::Subsection::Closed:: =#
    #= 8.2*Fresnel*integral*functions =#
    @apply_utils Antiderivative(FresnelS(~(a') + ~(b') * ~x), ~x) => ((~a + ~b * ~x) * FresnelS(~a + ~b * ~x)) / ~b + cos((pi / 2) * (~a + ~b * ~x) ^ 2) / (~b * pi) <-- FreeQ([~a, ~b], ~x)
    @apply_utils Antiderivative(FresnelC(~(a') + ~(b') * ~x), ~x) => ((~a + ~b * ~x) * FresnelC(~a + ~b * ~x)) / ~b - sin((pi / 2) * (~a + ~b * ~x) ^ 2) / (~b * pi) <-- FreeQ([~a, ~b], ~x)
    @apply_utils Antiderivative(FresnelS(~(a') + ~(b') * ~x) ^ 2, ~x) => ((~a + ~b * ~x) * FresnelS(~a + ~b * ~x) ^ 2) / ~b - 2 * Antiderivative((~a + ~b * ~x) * sin((pi / 2) * (~a + ~b * ~x) ^ 2) * FresnelS(~a + ~b * ~x), ~x) <-- FreeQ([~a, ~b], ~x)
    @apply_utils Antiderivative(FresnelC(~(a') + ~(b') * ~x) ^ 2, ~x) => ((~a + ~b * ~x) * FresnelC(~a + ~b * ~x) ^ 2) / ~b - 2 * Antiderivative((~a + ~b * ~x) * cos((pi / 2) * (~a + ~b * ~x) ^ 2) * FresnelC(~a + ~b * ~x), ~x) <-- FreeQ([~a, ~b], ~x)
    @apply_utils Antiderivative(FresnelS(~(a') + ~(b') * ~x) ^ ~n, ~x) => Unintegrable(FresnelS(~a + ~b * ~x) ^ ~n, ~x) <-- FreeQ([~a, ~b, ~n], ~x) && (NeQ(~n, 1) && NeQ(~n, 2))
    @apply_utils Antiderivative(FresnelC(~(a') + ~(b') * ~x) ^ ~n, ~x) => Unintegrable(FresnelC(~a + ~b * ~x) ^ ~n, ~x) <-- FreeQ([~a, ~b, ~n], ~x) && (NeQ(~n, 1) && NeQ(~n, 2))
    @apply_utils Antiderivative(FresnelS(~(b') * ~x) / ~x, ~x) => ((1 + im) / 4) * Antiderivative(erf((sqrt(pi) / 2) * (1 + im) * ~b * ~x) / ~x, ~x) + ((1 - im) / 4) * Antiderivative(erf((sqrt(pi) / 2) * (1 - im) * ~b * ~x) / ~x, ~x) <-- FreeQ(~b, ~x)
    @apply_utils Antiderivative(FresnelC(~(b') * ~x) / ~x, ~x) => ((1 - im) / 4) * Antiderivative(erf((sqrt(pi) / 2) * (1 + im) * ~b * ~x) / ~x, ~x) + ((1 + im) / 4) * Antiderivative(erf((sqrt(pi) / 2) * (1 - im) * ~b * ~x) / ~x, ~x) <-- FreeQ(~b, ~x)
    @apply_utils Antiderivative((~(d') * ~x) ^ ~(m') * FresnelS(~(b') * ~x), ~x) => ((~d * ~x) ^ (~m + 1) * FresnelS(~b * ~x)) / (~d * (~m + 1)) - (~b / (~d * (~m + 1))) * Antiderivative((~d * ~x) ^ (~m + 1) * sin((pi / 2) * (~b) ^ 2 * (~x) ^ 2), ~x) <-- FreeQ([~b, ~d, ~m], ~x) && NeQ(~m, -1)
    @apply_utils Antiderivative((~(d') * ~x) ^ ~(m') * FresnelC(~(b') * ~x), ~x) => ((~d * ~x) ^ (~m + 1) * FresnelC(~b * ~x)) / (~d * (~m + 1)) - (~b / (~d * (~m + 1))) * Antiderivative((~d * ~x) ^ (~m + 1) * cos((pi / 2) * (~b) ^ 2 * (~x) ^ 2), ~x) <-- FreeQ([~b, ~d, ~m], ~x) && NeQ(~m, -1)
    @apply_utils Antiderivative((~(c') + ~(d') * ~x) ^ ~(m') * FresnelS(~(a') + ~(b') * ~x), ~x) => ((~c + ~d * ~x) ^ (~m + 1) * FresnelS(~a + ~b * ~x)) / (~d * (~m + 1)) - (~b / (~d * (~m + 1))) * Antiderivative((~c + ~d * ~x) ^ (~m + 1) * sin((pi / 2) * (~a + ~b * ~x) ^ 2), ~x) <-- FreeQ([~a, ~b, ~c, ~d], ~x) && IGtQ(~m, 0)
    @apply_utils Antiderivative((~(c') + ~(d') * ~x) ^ ~(m') * FresnelC(~(a') + ~(b') * ~x), ~x) => ((~c + ~d * ~x) ^ (~m + 1) * FresnelC(~a + ~b * ~x)) / (~d * (~m + 1)) - (~b / (~d * (~m + 1))) * Antiderivative((~c + ~d * ~x) ^ (~m + 1) * cos((pi / 2) * (~a + ~b * ~x) ^ 2), ~x) <-- FreeQ([~a, ~b, ~c, ~d], ~x) && IGtQ(~m, 0)
    @apply_utils Antiderivative((~x) ^ ~(m') * FresnelS(~(b') * ~x) ^ 2, ~x) => ((~x) ^ (~m + 1) * FresnelS(~b * ~x) ^ 2) / (~m + 1) - ((2 * ~b) / (~m + 1)) * Antiderivative((~x) ^ (~m + 1) * sin((pi / 2) * (~b) ^ 2 * (~x) ^ 2) * FresnelS(~b * ~x), ~x) <-- FreeQ(~b, ~x) && (IntegerQ(~m) && NeQ(~m, -1))
    @apply_utils Antiderivative((~x) ^ ~(m') * FresnelC(~(b') * ~x) ^ 2, ~x) => ((~x) ^ (~m + 1) * FresnelC(~b * ~x) ^ 2) / (~m + 1) - ((2 * ~b) / (~m + 1)) * Antiderivative((~x) ^ (~m + 1) * cos((pi / 2) * (~b) ^ 2 * (~x) ^ 2) * FresnelC(~b * ~x), ~x) <-- FreeQ(~b, ~x) && (IntegerQ(~m) && NeQ(~m, -1))
    @apply_utils Antiderivative((~(c') + ~(d') * ~x) ^ ~(m') * FresnelS(~a + ~(b') * ~x) ^ 2, ~x) => (1 / (~b) ^ (~m + 1)) * Subst(Antiderivative(ExpandIntegrand(FresnelS(~x) ^ 2, ((~b * ~c - ~a * ~d) + ~d * ~x) ^ ~m, ~x), ~x), ~x, ~a + ~b * ~x) <-- FreeQ([~a, ~b, ~c, ~d], ~x) && IGtQ(~m, 0)
    @apply_utils Antiderivative((~(c') + ~(d') * ~x) ^ ~(m') * FresnelC(~a + ~(b') * ~x) ^ 2, ~x) => (1 / (~b) ^ (~m + 1)) * Subst(Antiderivative(ExpandIntegrand(FresnelC(~x) ^ 2, ((~b * ~c - ~a * ~d) + ~d * ~x) ^ ~m, ~x), ~x), ~x, ~a + ~b * ~x) <-- FreeQ([~a, ~b, ~c, ~d], ~x) && IGtQ(~m, 0)
    @apply_utils Antiderivative((~(c') + ~(d') * ~x) ^ ~(m') * FresnelS(~(a') + ~(b') * ~x) ^ ~(n'), ~x) => Unintegrable((~c + ~d * ~x) ^ ~m * FresnelS(~a + ~b * ~x) ^ ~n, ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~m, ~n], ~x)
    @apply_utils Antiderivative((~(c') + ~(d') * ~x) ^ ~(m') * FresnelC(~(a') + ~(b') * ~x) ^ ~(n'), ~x) => Unintegrable((~c + ~d * ~x) ^ ~m * FresnelC(~a + ~b * ~x) ^ ~n, ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~m, ~n], ~x)
    @apply_utils Antiderivative(ℯ ^ (~(c') + ~(d') * (~x) ^ 2) * FresnelS(~(b') * ~x), ~x) => ((1 + im) / 4) * Antiderivative(ℯ ^ (~c + ~d * (~x) ^ 2) * erf((sqrt(pi) / 2) * (1 + im) * ~b * ~x), ~x) + ((1 - im) / 4) * Antiderivative(ℯ ^ (~c + ~d * (~x) ^ 2) * erf((sqrt(pi) / 2) * (1 - im) * ~b * ~x), ~x) <-- FreeQ([~b, ~c, ~d], ~x) && EqQ((~d) ^ 2, (-(pi ^ 2) / 4) * (~b) ^ 4)
    @apply_utils Antiderivative(ℯ ^ (~(c') + ~(d') * (~x) ^ 2) * FresnelC(~(b') * ~x), ~x) => ((1 - im) / 4) * Antiderivative(ℯ ^ (~c + ~d * (~x) ^ 2) * erf((sqrt(pi) / 2) * (1 + im) * ~b * ~x), ~x) + ((1 + im) / 4) * Antiderivative(ℯ ^ (~c + ~d * (~x) ^ 2) * erf((sqrt(pi) / 2) * (1 - im) * ~b * ~x), ~x) <-- FreeQ([~b, ~c, ~d], ~x) && EqQ((~d) ^ 2, (-(pi ^ 2) / 4) * (~b) ^ 4)
    @apply_utils Antiderivative(ℯ ^ (~(c') + ~(d') * (~x) ^ 2) * FresnelS(~(a') + ~(b') * ~x) ^ ~(n'), ~x) => Unintegrable(ℯ ^ (~c + ~d * (~x) ^ 2) * FresnelS(~a + ~b * ~x) ^ ~n, ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~n], ~x)
    @apply_utils Antiderivative(ℯ ^ (~(c') + ~(d') * (~x) ^ 2) * FresnelC(~(a') + ~(b') * ~x) ^ ~(n'), ~x) => Unintegrable(ℯ ^ (~c + ~d * (~x) ^ 2) * FresnelC(~a + ~b * ~x) ^ ~n, ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~n], ~x)
    @apply_utils Antiderivative(sin(~(d') * (~x) ^ 2) * FresnelS(~(b') * ~x) ^ ~(n'), ~x) => ((pi * ~b) / (2 * ~d)) * Subst(Antiderivative((~x) ^ ~n, ~x), ~x, FresnelS(~b * ~x)) <-- FreeQ([~b, ~d, ~n], ~x) && EqQ((~d) ^ 2, (pi ^ 2 / 4) * (~b) ^ 4)
    @apply_utils Antiderivative(cos(~(d') * (~x) ^ 2) * FresnelC(~(b') * ~x) ^ ~(n'), ~x) => ((pi * ~b) / (2 * ~d)) * Subst(Antiderivative((~x) ^ ~n, ~x), ~x, FresnelC(~b * ~x)) <-- FreeQ([~b, ~d, ~n], ~x) && EqQ((~d) ^ 2, (pi ^ 2 / 4) * (~b) ^ 4)
    @apply_utils Antiderivative(sin(~c + ~(d') * (~x) ^ 2) * FresnelS(~(b') * ~x), ~x) => sin(~c) * Antiderivative(cos(~d * (~x) ^ 2) * FresnelS(~b * ~x), ~x) + cos(~c) * Antiderivative(sin(~d * (~x) ^ 2) * FresnelS(~b * ~x), ~x) <-- FreeQ([~b, ~c, ~d], ~x) && EqQ((~d) ^ 2, (pi ^ 2 / 4) * (~b) ^ 4)
    @apply_utils Antiderivative(cos(~c + ~(d') * (~x) ^ 2) * FresnelC(~(b') * ~x), ~x) => cos(~c) * Antiderivative(cos(~d * (~x) ^ 2) * FresnelC(~b * ~x), ~x) - sin(~c) * Antiderivative(sin(~d * (~x) ^ 2) * FresnelC(~b * ~x), ~x) <-- FreeQ([~b, ~c, ~d], ~x) && EqQ((~d) ^ 2, (pi ^ 2 / 4) * (~b) ^ 4)
    @apply_utils Antiderivative(sin(~(c') + ~(d') * (~x) ^ 2) * FresnelS(~(a') + ~(b') * ~x) ^ ~(n'), ~x) => Unintegrable(sin(~c + ~d * (~x) ^ 2) * FresnelS(~a + ~b * ~x) ^ ~n, ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~n], ~x)
    @apply_utils Antiderivative(cos(~(c') + ~(d') * (~x) ^ 2) * FresnelC(~(a') + ~(b') * ~x) ^ ~(n'), ~x) => Unintegrable(cos(~c + ~d * (~x) ^ 2) * FresnelC(~a + ~b * ~x) ^ ~n, ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~n], ~x)
    @apply_utils Antiderivative(cos(~(d') * (~x) ^ 2) * FresnelS(~(b') * ~x), ~x) => ((FresnelC(~b * ~x) * FresnelS(~b * ~x)) / (2 * ~b) - (1 / 8) * im * ~b * (~x) ^ 2 * HypergeometricPFQ([1, 1], [3 / 2, 2], (-1 / 2) * im * (~b) ^ 2 * pi * (~x) ^ 2)) + (1 / 8) * im * ~b * (~x) ^ 2 * HypergeometricPFQ([1, 1], [3 / 2, 2], (1 / 2) * im * (~b) ^ 2 * pi * (~x) ^ 2) <-- FreeQ([~b, ~d], ~x) && EqQ((~d) ^ 2, (pi ^ 2 / 4) * (~b) ^ 4)
    @apply_utils Antiderivative(sin(~(d') * (~x) ^ 2) * FresnelC(~(b') * ~x), ~x) => ((~b * pi * FresnelC(~b * ~x) * FresnelS(~b * ~x)) / (4 * ~d) + (1 / 8) * im * ~b * (~x) ^ 2 * HypergeometricPFQ([1, 1], [3 / 2, 2], -im * ~d * (~x) ^ 2)) - (1 / 8) * im * ~b * (~x) ^ 2 * HypergeometricPFQ([1, 1], [3 / 2, 2], im * ~d * (~x) ^ 2) <-- FreeQ([~b, ~d], ~x) && EqQ((~d) ^ 2, (pi ^ 2 / 4) * (~b) ^ 4)
    @apply_utils Antiderivative(cos(~c + ~(d') * (~x) ^ 2) * FresnelS(~(b') * ~x), ~x) => cos(~c) * Antiderivative(cos(~d * (~x) ^ 2) * FresnelS(~b * ~x), ~x) - sin(~c) * Antiderivative(sin(~d * (~x) ^ 2) * FresnelS(~b * ~x), ~x) <-- FreeQ([~b, ~c, ~d], ~x) && EqQ((~d) ^ 2, (pi ^ 2 / 4) * (~b) ^ 4)
    @apply_utils Antiderivative(sin(~c + ~(d') * (~x) ^ 2) * FresnelC(~(b') * ~x), ~x) => sin(~c) * Antiderivative(cos(~d * (~x) ^ 2) * FresnelC(~b * ~x), ~x) + cos(~c) * Antiderivative(sin(~d * (~x) ^ 2) * FresnelC(~b * ~x), ~x) <-- FreeQ([~b, ~c, ~d], ~x) && EqQ((~d) ^ 2, (pi ^ 2 / 4) * (~b) ^ 4)
    @apply_utils Antiderivative(cos(~(c') + ~(d') * (~x) ^ 2) * FresnelS(~(a') + ~(b') * ~x) ^ ~(n'), ~x) => Unintegrable(cos(~c + ~d * (~x) ^ 2) * FresnelS(~a + ~b * ~x) ^ ~n, ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~n], ~x)
    @apply_utils Antiderivative(sin(~(c') + ~(d') * (~x) ^ 2) * FresnelC(~(a') + ~(b') * ~x) ^ ~(n'), ~x) => Unintegrable(sin(~c + ~d * (~x) ^ 2) * FresnelC(~a + ~b * ~x) ^ ~n, ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~n], ~x)
    @apply_utils Antiderivative(~x * sin(~(d') * (~x) ^ 2) * FresnelS(~(b') * ~x), ~x) => (-(cos(~d * (~x) ^ 2)) * FresnelS(~b * ~x)) / (2 * ~d) + (1 / (2 * ~b * pi)) * Antiderivative(sin(2 * ~d * (~x) ^ 2), ~x) <-- FreeQ([~b, ~d], ~x) && EqQ((~d) ^ 2, (pi ^ 2 / 4) * (~b) ^ 4)
    @apply_utils Antiderivative(~x * cos(~(d') * (~x) ^ 2) * FresnelC(~(b') * ~x), ~x) => (sin(~d * (~x) ^ 2) * FresnelC(~b * ~x)) / (2 * ~d) - (~b / (4 * ~d)) * Antiderivative(sin(2 * ~d * (~x) ^ 2), ~x) <-- FreeQ([~b, ~d], ~x) && EqQ((~d) ^ 2, (pi ^ 2 / 4) * (~b) ^ 4)
    @apply_utils Antiderivative((~x) ^ ~m * sin(~(d') * (~x) ^ 2) * FresnelS(~(b') * ~x), ~x) => (-((~x) ^ (~m - 1)) * cos(~d * (~x) ^ 2) * FresnelS(~b * ~x)) / (2 * ~d) + (1 / (2 * ~b * pi)) * Antiderivative((~x) ^ (~m - 1) * sin(2 * ~d * (~x) ^ 2), ~x) + ((~m - 1) / (2 * ~d)) * Antiderivative((~x) ^ (~m - 2) * cos(~d * (~x) ^ 2) * FresnelS(~b * ~x), ~x) <-- FreeQ([~b, ~d], ~x) && (EqQ((~d) ^ 2, (pi ^ 2 / 4) * (~b) ^ 4) && IGtQ(~m, 1))
    @apply_utils Antiderivative((~x) ^ ~m * cos(~(d') * (~x) ^ 2) * FresnelC(~(b') * ~x), ~x) => (((~x) ^ (~m - 1) * sin(~d * (~x) ^ 2) * FresnelC(~b * ~x)) / (2 * ~d) - (~b / (4 * ~d)) * Antiderivative((~x) ^ (~m - 1) * sin(2 * ~d * (~x) ^ 2), ~x)) - ((~m - 1) / (2 * ~d)) * Antiderivative((~x) ^ (~m - 2) * sin(~d * (~x) ^ 2) * FresnelC(~b * ~x), ~x) <-- FreeQ([~b, ~d], ~x) && (EqQ((~d) ^ 2, (pi ^ 2 / 4) * (~b) ^ 4) && IGtQ(~m, 1))
    @apply_utils Antiderivative((~x) ^ ~m * sin(~(d') * (~x) ^ 2) * FresnelS(~(b') * ~x), ~x) => ((((~x) ^ (~m + 1) * sin(~d * (~x) ^ 2) * FresnelS(~b * ~x)) / (~m + 1) - (~d * (~x) ^ (~m + 2)) / (pi * ~b * (~m + 1) * (~m + 2))) + (~d / (pi * ~b * (~m + 1))) * Antiderivative((~x) ^ (~m + 1) * cos(2 * ~d * (~x) ^ 2), ~x)) - ((2 * ~d) / (~m + 1)) * Antiderivative((~x) ^ (~m + 2) * cos(~d * (~x) ^ 2) * FresnelS(~b * ~x), ~x) <-- FreeQ([~b, ~d], ~x) && (EqQ((~d) ^ 2, (pi ^ 2 / 4) * (~b) ^ 4) && ILtQ(~m, -2))
    @apply_utils Antiderivative((~x) ^ ~m * cos(~(d') * (~x) ^ 2) * FresnelC(~(b') * ~x), ~x) => ((((~x) ^ (~m + 1) * cos(~d * (~x) ^ 2) * FresnelC(~b * ~x)) / (~m + 1) - (~b * (~x) ^ (~m + 2)) / (2 * (~m + 1) * (~m + 2))) - (~b / (2 * (~m + 1))) * Antiderivative((~x) ^ (~m + 1) * cos(2 * ~d * (~x) ^ 2), ~x)) + ((2 * ~d) / (~m + 1)) * Antiderivative((~x) ^ (~m + 2) * sin(~d * (~x) ^ 2) * FresnelC(~b * ~x), ~x) <-- FreeQ([~b, ~d], ~x) && (EqQ((~d) ^ 2, (pi ^ 2 / 4) * (~b) ^ 4) && ILtQ(~m, -2))
    @apply_utils Antiderivative((~(e') * ~x) ^ ~(m') * sin(~(c') + ~(d') * (~x) ^ 2) * FresnelS(~(a') + ~(b') * ~x) ^ ~(n'), ~x) => Unintegrable((~e * ~x) ^ ~m * sin(~c + ~d * (~x) ^ 2) * FresnelS(~a + ~b * ~x) ^ ~n, ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~m, ~n], ~x)
    @apply_utils Antiderivative((~(e') * ~x) ^ ~(m') * cos(~(c') + ~(d') * (~x) ^ 2) * FresnelC(~(a') + ~(b') * ~x) ^ ~(n'), ~x) => Unintegrable((~e * ~x) ^ ~m * cos(~c + ~d * (~x) ^ 2) * FresnelC(~a + ~b * ~x) ^ ~n, ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~m, ~n], ~x)
    @apply_utils Antiderivative(~x * cos(~(d') * (~x) ^ 2) * FresnelS(~(b') * ~x), ~x) => (sin(~d * (~x) ^ 2) * FresnelS(~b * ~x)) / (2 * ~d) - (1 / (pi * ~b)) * Antiderivative(sin(~d * (~x) ^ 2) ^ 2, ~x) <-- FreeQ([~b, ~d], ~x) && EqQ((~d) ^ 2, (pi ^ 2 / 4) * (~b) ^ 4)
    @apply_utils Antiderivative(~x * sin(~(d') * (~x) ^ 2) * FresnelC(~(b') * ~x), ~x) => (-(cos(~d * (~x) ^ 2)) * FresnelC(~b * ~x)) / (2 * ~d) + (~b / (2 * ~d)) * Antiderivative(cos(~d * (~x) ^ 2) ^ 2, ~x) <-- FreeQ([~b, ~d], ~x) && EqQ((~d) ^ 2, (pi ^ 2 / 4) * (~b) ^ 4)
    @apply_utils Antiderivative((~x) ^ ~m * cos(~(d') * (~x) ^ 2) * FresnelS(~(b') * ~x), ~x) => (((~x) ^ (~m - 1) * sin(~d * (~x) ^ 2) * FresnelS(~b * ~x)) / (2 * ~d) - (1 / (pi * ~b)) * Antiderivative((~x) ^ (~m - 1) * sin(~d * (~x) ^ 2) ^ 2, ~x)) - ((~m - 1) / (2 * ~d)) * Antiderivative((~x) ^ (~m - 2) * sin(~d * (~x) ^ 2) * FresnelS(~b * ~x), ~x) <-- FreeQ([~b, ~d], ~x) && (EqQ((~d) ^ 2, (pi ^ 2 / 4) * (~b) ^ 4) && IGtQ(~m, 1))
    @apply_utils Antiderivative((~x) ^ ~m * sin(~(d') * (~x) ^ 2) * FresnelC(~(b') * ~x), ~x) => (-((~x) ^ (~m - 1)) * cos(~d * (~x) ^ 2) * FresnelC(~b * ~x)) / (2 * ~d) + (~b / (2 * ~d)) * Antiderivative((~x) ^ (~m - 1) * cos(~d * (~x) ^ 2) ^ 2, ~x) + ((~m - 1) / (2 * ~d)) * Antiderivative((~x) ^ (~m - 2) * cos(~d * (~x) ^ 2) * FresnelC(~b * ~x), ~x) <-- FreeQ([~b, ~d], ~x) && (EqQ((~d) ^ 2, (pi ^ 2 / 4) * (~b) ^ 4) && IGtQ(~m, 1))
    @apply_utils Antiderivative((~x) ^ ~m * cos(~(d') * (~x) ^ 2) * FresnelS(~(b') * ~x), ~x) => (((~x) ^ (~m + 1) * cos(~d * (~x) ^ 2) * FresnelS(~b * ~x)) / (~m + 1) - (~d / (pi * ~b * (~m + 1))) * Antiderivative((~x) ^ (~m + 1) * sin(2 * ~d * (~x) ^ 2), ~x)) + ((2 * ~d) / (~m + 1)) * Antiderivative((~x) ^ (~m + 2) * sin(~d * (~x) ^ 2) * FresnelS(~b * ~x), ~x) <-- FreeQ([~b, ~d], ~x) && (EqQ((~d) ^ 2, (pi ^ 2 / 4) * (~b) ^ 4) && ILtQ(~m, -1))
    @apply_utils Antiderivative((~x) ^ ~m * sin(~(d') * (~x) ^ 2) * FresnelC(~(b') * ~x), ~x) => (((~x) ^ (~m + 1) * sin(~d * (~x) ^ 2) * FresnelC(~b * ~x)) / (~m + 1) - (~b / (2 * (~m + 1))) * Antiderivative((~x) ^ (~m + 1) * sin(2 * ~d * (~x) ^ 2), ~x)) - ((2 * ~d) / (~m + 1)) * Antiderivative((~x) ^ (~m + 2) * cos(~d * (~x) ^ 2) * FresnelC(~b * ~x), ~x) <-- FreeQ([~b, ~d], ~x) && (EqQ((~d) ^ 2, (pi ^ 2 / 4) * (~b) ^ 4) && ILtQ(~m, -1))
    @apply_utils Antiderivative((~(e') * ~x) ^ ~(m') * cos(~(c') + ~(d') * (~x) ^ 2) * FresnelS(~(a') + ~(b') * ~x) ^ ~(n'), ~x) => Unintegrable((~e * ~x) ^ ~m * cos(~c + ~d * (~x) ^ 2) * FresnelS(~a + ~b * ~x) ^ ~n, ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~m, ~n], ~x)
    @apply_utils Antiderivative((~(e') * ~x) ^ ~(m') * sin(~(c') + ~(d') * (~x) ^ 2) * FresnelC(~(a') + ~(b') * ~x) ^ ~(n'), ~x) => Unintegrable((~e * ~x) ^ ~m * sin(~c + ~d * (~x) ^ 2) * FresnelC(~a + ~b * ~x) ^ ~n, ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~m, ~n], ~x)
    @apply_utils Antiderivative(FresnelS(~(d') * (~(a') + ~(b') * log(~(c') * (~x) ^ ~(n')))), ~x) => ~x * FresnelS(~d * (~a + ~b * log(~c * (~x) ^ ~n))) - ~b * ~d * ~n * Antiderivative(sin((pi / 2) * (~d * (~a + ~b * log(~c * (~x) ^ ~n))) ^ 2), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~n], ~x)
    @apply_utils Antiderivative(FresnelC(~(d') * (~(a') + ~(b') * log(~(c') * (~x) ^ ~(n')))), ~x) => ~x * FresnelC(~d * (~a + ~b * log(~c * (~x) ^ ~n))) - ~b * ~d * ~n * Antiderivative(cos((pi / 2) * (~d * (~a + ~b * log(~c * (~x) ^ ~n))) ^ 2), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~n], ~x)
    @apply_utils Antiderivative((~F)(~(d') * (~(a') + ~(b') * log(~(c') * (~x) ^ ~(n')))) / ~x, ~x) => (1 / ~n) * Subst((~F)(~d * (~a + ~b * ~x)), ~x, log(~c * (~x) ^ ~n)) <-- FreeQ([~a, ~b, ~c, ~d, ~n], ~x) && MemberQ([FresnelS, FresnelC], ~F)
    @apply_utils Antiderivative((~(e') * ~x) ^ ~(m') * FresnelS(~(d') * (~(a') + ~(b') * log(~(c') * (~x) ^ ~(n')))), ~x) => ((~e * ~x) ^ (~m + 1) * FresnelS(~d * (~a + ~b * log(~c * (~x) ^ ~n)))) / (~e * (~m + 1)) - ((~b * ~d * ~n) / (~m + 1)) * Antiderivative((~e * ~x) ^ ~m * sin((pi / 2) * (~d * (~a + ~b * log(~c * (~x) ^ ~n))) ^ 2), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~m, ~n], ~x) && NeQ(~m, -1)
    @apply_utils Antiderivative((~(e') * ~x) ^ ~(m') * FresnelC(~(d') * (~(a') + ~(b') * log(~(c') * (~x) ^ ~(n')))), ~x) => ((~e * ~x) ^ (~m + 1) * FresnelC(~d * (~a + ~b * log(~c * (~x) ^ ~n)))) / (~e * (~m + 1)) - ((~b * ~d * ~n) / (~m + 1)) * Antiderivative((~e * ~x) ^ ~m * cos((pi / 2) * (~d * (~a + ~b * log(~c * (~x) ^ ~n))) ^ 2), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~m, ~n], ~x) && NeQ(~m, -1)
end

max_stars_repo_path: src/rules/8 Special functions/8.2 Fresnel integral functions.jl


text:
using GeneralizedGenerated

export rand
EmptyNTtype = NamedTuple{(),Tuple{}} where T<:Tuple

rand(d::JointDistribution, N::Int) = [rand(d) for n in 1:N]

@inline function rand(m::JointDistribution)
    return _rand(getmoduletypencoding(m.model), m.model, m.args)
end

@inline function rand(m::Model)
    return _rand(getmoduletypencoding(m), m, NamedTuple())
end

@gg M function _rand(_::Type{M}, _m::Model, _args) where M <: TypeLevel{Module}
    Expr(:let,
        Expr(:(=), :M, from_type(M)),
        type2model(_m) |> sourceRand() |> loadvals(_args, NamedTuple()))
end

@gg M function _rand(_::Type{M}, _m::Model, _args::NamedTuple{()}) where M <: TypeLevel{Module}
    Expr(:let,
        Expr(:(=), :M, from_type(M)),
        type2model(_m) |> sourceRand())
end

export sourceRand
function sourceRand()
    function(m::Model)

        _m = canonical(m)
        proc(_m, st::Assign)  = :($(st.x) = $(st.rhs))
        proc(_m, st::Sample)  = :($(st.x) = rand($(st.rhs)))
        proc(_m, st::Return)  = :(return $(st.rhs))
        proc(_m, st::LineNumber) = nothing

        vals = map(x -> Expr(:(=), x,x),variables(_m))

        wrap(kernel) = @q begin
            $kernel
            $(Expr(:tuple, vals...))
        end

        buildSource(_m, proc, wrap) |> flatten
    end
end


text:
include("Graph.jl")
include("Planarity.jl")
include("VirtualEdge.jl")
include("PlanarDrawing.jl")

function DrawingAnalyze(G, Emb, X, Y, PATH; picSize = 20.0)
	ll = picSize / max(X[argmax(X)], Y[argmax(Y)])
	minAngle = 1e10

	vectorDot(p1, p2) = begin
		return p1[1]*p2[1] + p1[2]*p2[2]
	end
	vectorLen(p) = begin
		return sqrt(p[1]*p[1] + p[2]*p[2])
	end
	getNextIndex(x, y) = begin
		ret = (y == size(Emb[x], 1)) ? 1 : (y+1)
		while Emb[x][ret] > G.m
			ret = (ret == size(Emb[x], 1)) ? 1 : (ret+1)
		end
		return ret
	end
	getVector(u, eid) = begin
		if (X[u], Y[u]) == first(PATH[eid])
			return PATH[eid][2][1]-X[u], PATH[eid][2][2]-Y[u]
		end
		return PATH[eid][size(PATH[eid], 1)-1][1]-X[u], PATH[eid][size(PATH[eid], 1)-1][2]-Y[u]
	end
	getAngle(x, y) = begin
		if (y == getNextIndex(x, y))
			return 1e10
		end
		p1 = getVector(x, Emb[x][y])
		p2 = getVector(x, Emb[x][getNextIndex(x, y)])
		return acosd(vectorDot(p1, p2) / (vectorLen(p1) * vectorLen(p2)))
	end

	for i = 1 : G.n
		for j = 1 : size(Emb[i], 1)
			if Emb[i][j] <= G.m
				minAngle = min(minAngle, getAngle(i, j))
			end
		end
	end

	println("Minimum Angle : ", minAngle)

	getDistance(u, v) = begin
		return vectorLen((X[u]-X[v], Y[u]-Y[v]))
	end

	minP2Pdistance = 1e10
	for i = 1 : G.n, j = i+1 : G.n
		minP2Pdistance = min(minP2Pdistance, getDistance(i, j))
	end

	println("Minimum point-to-point distance : ", minP2Pdistance*ll)

	getiv(xa, ya, xb, yb, xx, yy) = begin
		lenab = vectorLen((xb-xa, yb-ya))
		ex = (xb - xa) / lenab
		ey = (yb - ya) / lenab
		vx = xx - xa
		vy = yy - ya
		pt = vx*ex + vy*ey
		xi = xa + ex*pt
		yi = ya + ey*pt
		return xi, yi
	end
	ivInAB(xi, yi, xa, ya, xb, yb) = begin
		if (xi < min(xa, xb)) || (xi > max(xa, xb))
			return false
		elseif (yi < min(ya, yb)) || (yi > max(ya, yb))
			return false
		end
		return true
	end
	getP2Ldistance(p, la, lb) = begin
		xi, yi = getiv(la[1], la[2], lb[1], lb[2], p[1], p[2])
		if !ivInAB(xi, yi, la[1], la[2], lb[1], lb[2])
			return 1e10
		end
		return vectorLen((p[1]-xi, p[2]-yi))
	end

	minP2Ldistance = 1e10

	updateP2L(x) = begin
		if (x < minP2Ldistance)
			minP2Ldistance = x
		end
	end

	for i = 1 : G.n, j = 1 : G.m
		if (G.E[j][2] != i) && (G.E[j][3] != i)
			for k = 1 : size(PATH[j], 1)-1
				updateP2L(getP2Ldistance((X[i], Y[i]), PATH[j][k], PATH[j][k+1]))
			end
		end
	end

	println("Minimum point-to-line distance : ", minP2Ldistance*ll)
end

function PlanarPolyDraw(G)
	# get planar subgraph
	edgeIndex = zeros(Int, G.m)
	foreach(i -> edgeIndex[i] = i, 1 : G.m)
	sort!(edgeIndex, by = x->-G.E[x][4])
	EdgeList = Array{Tuple{Int32, Int32, Int32, Float64}, 1}()
	preservedEdgeList = Array{Tuple{Int32, Int32}, 1}()
	deletedEdgeList = Array{Tuple{Int32, Int32}, 1}()
	deletedEdgeC = Array{Float64, 1}()
	preservedEdgeC = 0.0
	M2 = 0
	for i = 1 : G.m
		(ID, u, v, w) = G.E[edgeIndex[i]]
		M2 += 1
		push!(EdgeList, (M2, u, v, 1.0))
		if LRtest(Graph(G.n, size(EdgeList, 1), G.V, EdgeList)) == nothing
			M2 -= 1
			pop!(EdgeList)
			push!(deletedEdgeList, (u, v))
			push!(deletedEdgeC, w)
		else
			preservedEdgeC += w
			push!(preservedEdgeList, (u, v))
		end
	end

	# add some virtual edge to make the graph connected
	EdgeList = copy(AddToConnected(Graph(G.n, size(EdgeList, 1), G.V, EdgeList)))
	# add some virtual edge to make the graph biconnected
	Emb = copy(LRtest(Graph(G.n, size(EdgeList, 1), G.V, EdgeList), needEmbedding = true))
	EdgeList = copy(AddToBiConnected(Graph(G.n, size(EdgeList, 1), G.V, EdgeList), Emb))
	# add some virtual edge to get a maximal planar graph
	Emb = copy(LRtest(Graph(G.n, size(EdgeList, 1), G.V, EdgeList), needEmbedding = true))
	EdgeList = copy(AddToMaximalPlanar(Graph(G.n, size(EdgeList, 1), G.V, EdgeList), Emb))
	# calculate canonical ordering
	Emb = copy(LRtest(Graph(G.n, size(EdgeList, 1), G.V, EdgeList), needEmbedding = true))
	Order = CanonicalOrdering(Graph(G.n, size(EdgeList, 1), G.V, EdgeList), Emb)

	# using polyline drawing to draw a planar graph
	X, Y, PATH, L, H = PolylineDrawing(Graph(G.n, size(EdgeList, 1), G.V, EdgeList), M2, Order, Emb)

	#DrawingAnalyze(Graph(G.n, M2, G.V, EdgeList[1:M2]), Emb, X, Y, PATH)

	return X, Y, PATH, deletedEdgeList, L, H, preservedEdgeC, deletedEdgeC
end


text:
using Base.Broadcast
using Base.Broadcast: broadcasted, Broadcasted

# Forward mode broadcast rule
struct FwdBroadcast{N, T<:AbstractTangentBundle{N}}
  f::T
end
(f::FwdBroadcast{N})(args::AbstractTangentBundle{N}...) where {N} = ∂☆{N}()(f.f, args...)

n_getfield(∂ₙ::∂☆{N}, b::ATB{N}, x::Union{Symbol, Int}) where {N} = ∂ₙ(ZeroBundle{N}(getfield), b, ZeroBundle{N}(x))

function (∂ₙ::∂☆{N})(zc::ZeroBundle{N, typeof(copy)},
                     bc::ATB{N, <:Broadcasted}) where {N}
  bc = ∂ₙ(ZeroBundle{N}(Broadcast.flatten), bc)
  args = n_getfield(∂ₙ, bc, :args)
  r = copy(Broadcasted(
      FwdMap(n_getfield(∂ₙ, bc, :f)),
      ntuple(length(primal(args))) do i
        val = n_getfield(∂ₙ, args, i)
        if ndims(primal(val)) == 0
          return Ref(∂ₙ(ZeroBundle{N}(getindex), val))
        else
          return unbundle(val)
        end
      end))
  if isa(r, AbstractArray)
    r = rebundle(r)
  end
  return r
end

# Broadcast over one element is just map
function (∂⃖ₙ::∂⃖{N})(::typeof(broadcasted), f, a::Array) where {N}
    ∂⃖ₙ(map, f, a)
end

# The below is from Zygote: TODO: DO we want to do something better here?

accum_sum(xs::Nothing; dims = :) = NoTangent()
accum_sum(xs::AbstractArray{Nothing}; dims = :) = NoTangent()
accum_sum(xs::AbstractArray{<:Number}; dims = :) = sum(xs, dims = dims)
accum_sum(xs::AbstractArray{<:AbstractArray{<:Number}}; dims = :) = sum(xs, dims = dims)
accum_sum(xs::Number; dims = :) = xs

# https://github.com/FluxML/Zygote.jl/issues/594
function Base.reducedim_init(::typeof(identity), ::typeof(accum), A::AbstractArray, region)
  Base.reducedim_initarray(A, region, NoTangent(), Union{Nothing,eltype(A)})
end

trim(x, Δ) = reshape(Δ, ntuple(i -> size(Δ, i), Val(ndims(x))))

unbroadcast(x::AbstractArray, x̄) =
  size(x) == size(x̄) ? x̄ :
  length(x) == length(x̄) ? trim(x, x̄) :
    trim(x, accum_sum(x̄, dims = ntuple(i -> size(x, i) == 1 ? i : ndims(x̄)+1, Val(ndims(x̄)))))

unbroadcast(x::Number, x̄) = accum_sum(x̄)
unbroadcast(x::Tuple{<:Any}, x̄) = (accum_sum(x̄),)
unbroadcast(x::Base.RefValue, x̄) = (x=accum_sum(x̄),)

unbroadcast(x::AbstractArray, x̄::Nothing) = NoTangent()

const Numeric = Union{Number, AbstractArray{<:Number, N} where N}

function ChainRulesCore.rrule(::typeof(broadcasted), ::typeof(+), xs::Numeric...)
    broadcast(+, xs...), ȳ -> (NoTangent(), NoTangent(), map(x -> unbroadcast(x, unthunk(ȳ)), xs)...)
end

ChainRulesCore.rrule(::typeof(broadcasted), ::typeof(-), x::Numeric, y::Numeric) = x .- y,
  Δ -> let Δ=unthunk(Δ); (NoTangent(), NoTangent(), unbroadcast(x, Δ), -unbroadcast(y, Δ)); end

ChainRulesCore.rrule(::typeof(broadcasted), ::typeof(*), x::Numeric, y::Numeric) = x.*y,
  z̄ -> let z̄=unthunk(z̄); (NoTangent(), NoTangent(), unbroadcast(x, z̄ .* conj.(y)), unbroadcast(y, z̄ .* conj.(x))); end