text:
abstract type AbstractCamera end

mutable struct FreeCamera <: AbstractCamera
	deltaTrans::Vector{Float32}
	deltaAngles::Vector{Float32}

	FreeCamera(deltaTrans::Vector{Float32}, deltaAngles::Vector{Float32}) = new(deltaTrans, deltaAngles)
end

function update(cam::GRU.Camera, trans::Vector{Float32}, angles::Vector{Float32})
	m = Matrix{Float32}(I, 4, 4)
	Math3D.rotxyz(m, angles...)
	Math3D.trans(m, trans)
	GRU.settransform(cam, GRU.gettransform(cam)*m)
end

function process_input(engine::Engine, cam::FreeCamera)
	trans = Array{Float32}(undef, 3)
	angles = Array{Float32}(undef, 3)
	deltaTime = deltatime(engine)

	trans[1] = (GLFW.GetKey(engine.window, GLFW.KEY_D) - GLFW.GetKey(engine.window, GLFW.KEY_A)) * cam.deltaTrans[1] * deltaTime
	trans[2] = (GLFW.GetKey(engine.window, GLFW.KEY_F) - GLFW.GetKey(engine.window, GLFW.KEY_R)) * cam.deltaTrans[2] * deltaTime
	trans[3] = (GLFW.GetKey(engine.window, GLFW.KEY_W) - GLFW.GetKey(engine.window, GLFW.KEY_S)) * cam.deltaTrans[3] * deltaTime

	angles[1] = (GLFW.GetKey(engine.window, GLFW.KEY_G) - GLFW.GetKey(engine.window, GLFW.KEY_T)) * cam.deltaAngles[1] * deltaTime
	angles[2] = (GLFW.GetKey(engine.window, GLFW.KEY_E) - GLFW.GetKey(engine.window, GLFW.KEY_Q)) * cam.deltaAngles[2] * deltaTime
	angles[3] = (GLFW.GetKey(engine.window, GLFW.KEY_Z) - GLFW.GetKey(engine.window, GLFW.KEY_C)) * cam.deltaAngles[3] * deltaTime

	update(engine.renderer.camera, trans, angles)
end


text:
module AdvancedHMC

const DEBUG = Bool(parse(Int, get(ENV, "DEBUG_AHMC", "0")))

using Statistics: mean, var, middle
using LinearAlgebra: Symmetric, UpperTriangular, mul!, ldiv!, dot, I, diag, cholesky
using LazyArrays: BroadcastArray
using Random: GLOBAL_RNG, AbstractRNG

import StatsBase: sample

# Notations
# d - dimension of sampling sapce
# π(θ) - target distribution
# r - momentum variable

include("metric.jl")
export UnitEuclideanMetric, DiagEuclideanMetric, DenseEuclideanMetric
include("hamiltonian.jl")
export Hamiltonian
include("integrator.jl")
export Leapfrog
include("proposal.jl")
export TakeLastProposal, find_good_eps, NUTS

include("adaptation.jl")
export NesterovDualAveraging, PreConditioner, NaiveCompAdaptor, StanNUTSAdaptor
include("diagnosis.jl")
include("sampler.jl")
export sample

end # module

ext: jl


text:
using CRC32
using Base.Test

@test crc32("") == 0
@test crc32("hello") == 0x3610a686
@test crc32("Julia programming language") == 0xfc485364

crc = crc32("programming")
@test crc == 0x41f292cf
@test crc32("programming language") == 0xe83fb5fd
@test crc32(" language", crc) == 0xe83fb5fd


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

__precompile__(true)

"""
    MGVI

An implementation of the Metric Gaussian Variational Inference algorithm.
"""
module MGVI

import ChainRulesCore
using Distributed
using LinearAlgebra
using Random
using SparseArrays
using Base.Iterators
using Distributions
using DistributionsAD
using FillArrays
import ForwardDiff
using LineSearches
using LinearMaps
using IterativeSolvers
using Optim
using Parameters
using PDMats
using PositiveFactorizations
import SparseArrays: blockdiag
using SparseArrays
using StaticArrays
using ValueShapes
import Zygote
using DocStringExtensions

using ChainRulesCore: AbstractTangent, Tangent, NoTangent, ZeroTangent, ProjectTo, AbstractThunk, unthunk
import Statistics: mean

include("util.jl")
include("custom_linear_maps.jl")
include("shapes.jl")
include("jacobian_maps.jl")
include("information.jl")
include("residual_samplers.jl")
include("newtoncg.jl")
include("mgvi_impl.jl")

end # module


text:
#Worked with Kasra and Waleed#

cd("D:\\Semester 3_Fall 2020\\ECON 6453 - Econometrics III\\ProblemSets\\PS2-optimization-intro")
using Random, Distributions, JLD2, DataFrames, CSV, LinearAlgebra, Statistics, FreqTables, TexTables, Optim, HTTP, GLM

#:::::::::::::::::::::::::::::::::::::::::::::::::::
# question 6: Wrapping Around a Function#
#:::::::::::::::::::::::::::::::::::::::::::::::::::
function ProblemSet_2()
#:::::::::::::::::::::::::::::::::::::::::::::::::::
# question 1
#:::::::::::::::::::::::::::::::::::::::::::::::::::
f(x) = -x[1]^4-10x[1]^3-2x[1]^2-3x[1]-2
minusf(x) = x[1]^4+10x[1]^3+2x[1]^2+3x[1]+2
startval = rand(1)   # random starting value
result = optimize(minusf, startval, LBFGS())

#:::::::::::::::::::::::::::::::::::::::::::::::::::
# question 2
#:::::::::::::::::::::::::::::::::::::::::::::::::::
url = "https://raw.githubusercontent.com/OU-PhD-Econometrics/fall-2020/master/ProblemSets/PS1-julia-intro/nlsw88.csv"
df = CSV.read(HTTP.get(url).body)
X = [ones(size(df,1),1) df.age df.race.==1 df.collgrad.==1]
y = df.married.==1

function ols(beta, X, y)
    ssr = (y.-X*beta)'*(y.-X*beta)
    return ssr
end

beta_hat_ols = optimize(b -> ols(b, X, y), rand(size(X,2)), LBFGS(), Optim.Options(g_tol=1e-6, iterations=100_000, show_trace=true))
println(beta_hat_ols.minimizer)

bols = inv(X'*X)*X'*y
df.white = df.race.==1
bols_lm = lm(@formula(married ~ age + white + collgrad), df)

#:::::::::::::::::::::::::::::::::::::::::::::::::::
# question 3: Using Optim for Logit Regression#
#:::::::::::::::::::::::::::::::::::::::::::::::::::
function logit(alpha, X, d)
    loglike_function=[d[i]*log((1+exp((X*alpha)[i])))+(1-d[i])*log(1-(1+exp((X*alpha)[i]))) for i in 1:size(d,1)]
    # P_[i1]=1+exp(X*alpha) #
    loglike=-sum(loglike_function)
    return loglike
end

alpha_loglike = optimize(b -> logit(b, X, y), rand(size(X,2)), LBFGS(), Optim.Options(g_tol=1e-6, iterations=100_000, show_trace=true))
println(alpha_loglike.minimizer)

#:::::::::::::::::::::::::::::::::::::::::::::::::::
# question 4: Using GLM for Logit Regression#
#:::::::::::::::::::::::::::::::::::::::::::::::::::
df.white = df.race.==1
α_glm = glm(@formula(married ~ age + white + collgrad), df, Binomial(), LogitLink())
println(α_glm)

#:::::::::::::::::::::::::::::::::::::::::::::::::::
# question 5: Multinomial Logistic Regression#
#:::::::::::::::::::::::::::::::::::::::::::::::::::
freqtable(df, :occupation) # note small number of obs in some occupations
df = dropmissing(df, :occupation)
df[df.occupation.==10,:occupation] .= 9
df[df.occupation.==11,:occupation] .= 9
df[df.occupation.==12,:occupation] .= 9
df[df.occupation.==13,:occupation] .= 9

X = [ones(size(df,1),1) df.age df.race.==1 df.collgrad.==1]
y = df.occupation

function mnlogit(alpha, x, d)
    d=[i==j ? 1 : 0 for i in d , j in 1:6]
    loglike_function2=[log(exp.(x[i,:]'*alpha[:,j])/(1+sum(exp.(x[i,:]'*alpha)))) for i in 1:size(x,1) ,j in 1:6]
    loglike_multi=-sum(d.*loglike_function2)
    return loglike_multi
end

alpha_mnlog = optimize(b -> mnlogit(b, X, y), rand(size(X,2),6), LBFGS(), Optim.Options(g_tol=1e-7, iterations=100_000, show_trace=true))
println(alpha_mnlog.minimizer)

#:::::::::::::::::::::::::::::::::::::::::::::::::::
# question 6: Wrapping Around a Function#
#:::::::::::::::::::::::::::::::::::::::::::::::::::
end
ProblemSet_2()

                            #The End#


text:
"""
    SeisEnvelope(in; <keyword arguments>)

Calculate the envelope attribute of a group of input traces.

# Arguments
- `in`: input data. A 2D Array where the first dimension is time.

# Example
```julia
julia> using PyPlot
julia> dtsec = 0.002; w = Ricker(dt=dtsec);
julia> e = SeisEnvelope(w); t=dtsec*collect(0:1:length(w)-1);plot(t,w,t,e,"-r");xlabel("Time [s]")
```
*Credits: Aaron Stanton,2017*

"""
function SeisEnvelope(d)

	D = fft(d,1)
	D[1:Int(floor(size(d,1)/2)),:] .= 0.0
	return 2*abs.(ifft(D,1))

end


text:
using LinearAlgebra

include("../quadrature/quadgrid.jl")

function param_grid_comp(n_param, param_priorInfo)
    # extract param info
    thetamin = param_priorInfo.theta.min
    thetamax = param_priorInfo.theta.max
    lambdamin = param_priorInfo.lambda.min
    lambdamax = param_priorInfo.lambda.max
    n_theta = n_param.theta
    n_lambda = n_param.lambda

   # form a 2d theta grid
   theta1_grid, wt1_theta = quadgrid(thetamin[1], thetamax[1], n_theta[1])
   theta2_grid, wt2_theta = quadgrid(thetamin[2], thetamax[2], n_theta[2])
   n_theta_totl = n_theta[1] * n_theta[2]
   theta_grid = zeros(n_theta_totl, 2)
   wt_theta = zeros(n_theta_totl)
   for i in 1:n_theta[1], j in 1:n_theta[2]
       theta_grid[j + (i-1) * n_theta[2], :] = [theta1_grid[i] theta2_grid[j]]
       wt_theta[j + (i-1) * n_theta[2]] = wt1_theta[i] * wt2_theta[j]
   end

   # form a 1d lambda grid
   lambda_grid, wt_lambda = quadgrid(lambdamin, lambdamax, n_lambda)

   thetaInfo = (grid = theta_grid, wt = wt_theta)
   lambdaInfo = (grid = lambda_grid, wt = wt_lambda)
   param_gridInfo = (theta = thetaInfo, lambda = lambdaInfo)

   return param_gridInfo
end

max_stars_repo_name: dbindel/btg

max_issues_repo_name: dbindel/btg

max_forks_repo_name: dbindel/btg


text:
using Test

using DualOptimalFiltering
# Run tests
@time @test 1 == 1

include("test_common_utility_functions.jl")
include("test_mcmc_sampler.jl")
include("test_pruning_functions.jl")
include("test_kde_for_pf_samples.jl")
include("test_dirichlet_kde.jl")
include("test_exact_L2_distances.jl")
include("test_post_process_Dirichlet_mixture_posterior.jl")

include("test_functions_CIR.jl")
include("test_filtering_CIR.jl")
include("test_approx_filtering_CIR.jl")
include("test_likelihood_CIR.jl")
include("test_smoothing_CIR.jl")
include("test_joint_smoothing_CIR.jl")
include("test_approximate_smoothing_CIR.jl")
include("test_reparam_CIR.jl")
include("test_full_inference_CIR.jl")

include("test_approx_filtering_WF.jl")
include("test_likelihood_WF.jl")
include("test_particle_filtering_WF.jl")

max_stars_repo_name: konkam/DualOptimalFiltering

max_issues_repo_name: konkam/DualOptimalFiltering

max_forks_repo_name: konkam/DualOptimalFiltering


text:
# https://wickedengine.net/2020/04/26/capsule-collision-detection/

struct CapsuleContact{T}
    position::Vector{T}
    orientation::Quaternion{T}
    height::T
    radius::T
end

p1 = [0.0, 0.0, 0.0]
# o1 = one(Quaternion)
o1 = Quaternion(RotX(0.6) * RotY(0.0 * π))
h1 = 1.0
r1 = 0.1
cap1 = CapsuleContact(p1, o1, h1, r1)

p2 = [1.0, 0.0, 0.0]
o2 = one(Quaternion)
o2 = Quaternion(RotZ(0.0) * RotY(-0.25 * π))
h2 = 1.0
r2 = 0.1
cap2 = CapsuleContact(p2, o2, h2, r2)

function get_tip_base(cap::CapsuleContact)
    k = [0.0, 0.0, 0.5 * cap.height + cap.radius]
    tip = cap.position + vector_rotate(k, cap.orientation)
    base = cap.position + vector_rotate(-k, cap.orientation)
    return Array(tip), Array(base)
end

t1, b1 = get_tip_base(cap1)
t2, b2 = get_tip_base(cap2)

function closest_point_on_line(pa::Vector{T}, pb::Vector{T}, p::Vector{T}) where T
    pab = pb - pa
    t = dot(p - pa, pab) / dot(pab, pab)
    return pa + min(max(t, 0.0), 1.0) * pab
end

function closest_segment_points_capsule_capsule(capa::CapsuleContact, capb::CapsuleContact)
    # capsule a
    a_tip, a_base = get_tip_base(capa)
    a_radius = capa.radius

    # capsule b
    b_tip, b_base = get_tip_base(capb)
    b_radius = capb.radius

    # capsule a
    a_normal = a_tip - a_base
    a_normal ./= norm(a_normal)
    a_line_end_offset = a_normal * a_radius
    a_A = a_base + a_line_end_offset
    a_B = a_tip - a_line_end_offset

    # capsule b
    b_normal = b_tip - b_base
    b_normal ./= norm(b_normal)
    b_line_end_offset = b_normal * b_radius
    b_A = b_base + b_line_end_offset
    b_B = b_tip - b_line_end_offset

    # vectors between line endpoints
    v0 = b_A - a_A
    v1 = b_B - a_A
    v2 = b_A - a_B
    v3 = b_B - a_B

    # squared distances
    d0 = dot(v0, v0)
    d1 = dot(v1, v1)
    d2 = dot(v2, v2)
    d3 = dot(v3, v3)

    # best potential endpoint on capsule a
    if (d2 < d0 || d2 < d1 || d3 < d0 || d3 < d1)
        bestA = a_B
    else
        bestA = a_A
    end

    # select point on capsule B line segment nearest best potential endpoint on A capsule
    bestB = closest_point_on_line(b_A, b_B, bestA)

    # same for capsule A segment
    bestA = closest_point_on_line(a_A, a_B, bestB)

    return bestA, bestB
    # # direction
    # direction = bestB - bestA
    # direction_magnitude = norm(direction)
    # direction ./= direction_magnitude

    # c1 = bestA + direction * a_radius
    # c2 = bestB - direction * b_radius
    # # distance = direction_magnitude - (a_radius + b_radius)
    # return c1, c2
end

function kinematics_capsule_capsule(capa::CapsuleContact, capb::CapsuleContact)
    bestA, bestB = closest_segment_points_capsule_capsule(capa, capb)

    kA = vector_rotate(bestA - capa.position, inv(capa.orientation))
    kB = vector_rotate(bestB - capb.position, inv(capb.orientation))

    return kA, kB
end

ka, kb = kinematics_capsule_capsule(cap1, cap2)

function closest_points_capsule_capsule(capa::CapsuleContact, capb::CapsuleContact)
    @show "hi"
    bestA, bestB = closest_segment_points_capsule_capsule(capa, capb)

    direction = bestB - bestA
    direction_magnitude = norm(direction)
    direction ./= direction_magnitude

    c1 = bestA + direction * capa.radius
    c2 = bestB - direction * capb.radius
    # distance = direction_magnitude - (a_radius + b_radius)
    return c1, c2
end

function direction_capsule_capsule(capa::CapsuleContact, capb::CapsuleContact)
    ca, cb = closest_segment_points_capsule_capsule(capa, capb)

    direction = cb - ca
    direction_magnitude = norm(direction)

    if direction_magnitude > 0.0
        direction ./= direction_magnitude
    else
        direction = [1.0, 0.0, 0.0]
    end

    return direction
end

function sdf_capsule_capsule(capa::CapsuleContact, capb::CapsuleContact)
    ca, cb = closest_points_capsule_capsule(capa, capb)

    direction = cb - ca
    direction_magnitude = norm(direction)

    return direction_magnitude
end

c1, c2 = closest_points_capsule_capsule(cap1, cap2)

# dir = direction_capsule_capsule(cap1, cap2)

# dist = distance_capsule_capsule(cap1, cap2)

# pa = [0.0, 0.0, 0.0]
# pb = [1.0, 1.0, 0.0]
# p = [1.0, 0.5, 0.0]

# closest_point_on_line(pa, pb, p)

# # capsule a
# a_tip = [1.0, 1.0, 0.0]
# a_base = [0.0, 0.0, 0.0]
# a_radius = 0.01

# # capsule b
# b_tip = [2.0, 1.0, 0.0]
# b_base = [2.0, 0.0, 0.0]
# b_radius = 0.01

# # capsule a
# a_normal = a_tip - a_base
# a_normal ./= norm(a_normal)
# a_line_end_offset = a_normal * a_radius
# a_A = a_base + a_line_end_offset
# a_B = a_tip - a_line_end_offset

# # capsule b
# b_normal = b_tip - b_base
# b_normal ./= norm(b_normal)
# b_line_end_offset = b_normal * b_radius
# b_A = b_base + b_line_end_offset
# b_B = b_tip - b_line_end_offset

# # vectors between line endpoints
# v0 = b_A - a_A
# v1 = b_B - a_A
# v2 = b_A - a_B
# v3 = b_B - a_B

# # squared distances
# d0 = dot(v0, v0)
# d1 = dot(v1, v1)
# d2 = dot(v2, v2)
# d3 = dot(v3, v3)

# # best potential endpoint on capsule a
# if (d2 < d0 || d2 < d1 || d3 < d0 || d3 < d1)
#     bestA = a_B
# else
#     bestA = a_A
# end
# # select point on capsule B line segment nearest best potential endpoint on A capsule
# bestB = closest_point_on_line(b_A, b_B, bestA)

# # same for capsule A segment
# bestA = closest_point_on_line(a_A, a_B, bestB)

# # direction
# direction = bestB - bestA
# direction_magnitude = norm(direction)
# direction ./= direction_magnitude

# distance = direction_magnitude - (a_radius + b_radius)

# visualize
vis = Visualizer()
open(vis)

# capsule a
cyl1 = GeometryBasics.Cylinder(Point(0.0, 0.0, -0.5 * cap1.height), Point(0.0, 0.0, 0.5 * cap1.height), cap1.radius)
sph1 = GeometryBasics.Sphere(Point(0.0, 0.0, 0.0), cap1.radius)
color1 = Colors.RGBA(0.7, 0.7, 0.7, 0.5);
setobject!(vis[:cap1][:cyl], cyl1, MeshPhongMaterial(color=color1))
setobject!(vis[:cap1][:base], sph1, MeshPhongMaterial(color=color1))
settransform!(vis[:cap1][:base], Translation(0.0, 0.0, -0.5 * cap1.height))
setobject!(vis[:cap1][:tip], sph1, MeshPhongMaterial(color=color1))
settransform!(vis[:cap1][:tip], Translation(0.0, 0.0, 0.5 * cap1.height))

# capsule b
cyl2 = GeometryBasics.Cylinder(Point(0.0, 0.0, -0.5 * cap2.height), Point(0.0, 0.0, 0.5 * cap2.height), cap2.radius)
sph2 = GeometryBasics.Sphere(Point(0.0, 0.0, 0.0), cap2.radius)
color2 = Colors.RGBA(0.7, 0.7, 0.7, 0.5);
setobject!(vis[:cap2][:cyl], cyl2, MeshPhongMaterial(color=color2))
setobject!(vis[:cap2][:base], sph2, MeshPhongMaterial(color=color2))
settransform!(vis[:cap2][:base], Translation(0.0, 0.0, -0.5 * cap2.height))
setobject!(vis[:cap2][:tip], sph2, MeshPhongMaterial(color=color2))
settransform!(vis[:cap2][:tip], Translation(0.0, 0.0, 0.5 * cap2.height))

# set configuration
settransform!(vis[:cap1], compose(Translation(cap1.position), LinearMap(cap1.orientation)))
settransform!(vis[:cap2], compose(Translation(cap2.position), LinearMap(cap2.orientation)))

# closest points
color_cp = Colors.RGBA(0.0, 0.0, 0.0, 1.0);
cs1 = GeometryBasics.Sphere(Point(0.0, 0.0, 0.0), 0.025)
cs2 = GeometryBasics.Sphere(Point(0.0, 0.0, 0.0), 0.025)
setobject!(vis[:cs1], cs1, MeshPhongMaterial(color=color_cp))
setobject!(vis[:cs2], cs2, MeshPhongMaterial(color=color_cp))
settransform!(vis[:cs1], Translation(c1))
settransform!(vis[:cs2], Translation(c2))


text:
module Analog

export find_point, error_cov

using Statistics

using NearestNeighbors

function find_point(r, tree, p, k, f)
    ind, dist = knn(tree, p[:], k)
    mask = (ind .+ f) .<= size(r)[1]
    dist = 1 ./ dist[mask]
    ind = ind[mask]
    return sum(dist .* r[ind .+ f, :], dims=1)/sum(dist)
end

function error_cov(y, r, M, window, k, k_r, osc_vars; validation_pct=0.1)
    validation = round(Int, validation_pct*size(y)[1])
    tree = KDTree(copy((y[validation+1:end, :])'))
    tree_r = KDTree(copy((r[validation+1:end, :])'))
    errs = Array{Float64}(undef, length(osc_vars), length(M:validation-window))

    for (i, i_p) in enumerate(M:validation-window)
        p = y[i_p, :]
        p_r = find_point(r, tree, p, k, validation)

        forecast = find_point(r, tree_r, p_r[:], k_r, validation + window)
        err = r[i_p + window, :] - forecast'

        errs[:, i] = err
    end

    R = cov(errs')
    return R
end

function error_corr(y, r, M, window, k, k_r, osc_vars; validation_pct=0.1)
    validation = round(Int, validation_pct*size(y)[1])
    tree = KDTree(copy((y[validation+1:end, osc_vars])'))
    tree_r = KDTree(copy((r[validation+1:end, :])'))
    corrs = Array{Float64}(undef, length(M:validation-window))
    errs_m = Array{Float64}(undef, length(M:validation-window))
    m = mean(r, dims=1)

    for (i, i_p) in enumerate(M:validation-window)
        p = y[i_p, osc_vars]
        p_r = find_point(r, tree, p, k, validation)

        forecast = find_point(r, tree_r, p_r[:], k_r, validation + window)

        corrs[i] = sqrt(mean((r[i_p + window, :] - (forecast')[:]).^2))
        errs_m[i] = sqrt(mean((r[i_p + window, :] - m').^2))
    end

    return mean(corrs), mean(errs_m)
end

end


text:
using ResumableFunctions

extension(url::String) = try
    match(r"\.[A-Za-z0-9]+$", url).match
catch
    ""
end

@resumable function walkfiles(root; exts=Set((".md", ".html")),
                              ex_dirs::AbstractSet=Set(),
                              dirs::Union{AbstractSet,Nothing}=nothing,
                              subdir=false)
    """
iterate over files in a directory, recursively and selectively by extension name and dir name
"""
    for p in readdir(root)
        path = joinpath(root, p)
        # directory, not excluded
        if isdir(path)
            name = splitpath(p)[end]
            if in(name, ex_dirs)
                continue # requires ResumableFunctions > 0.0.6
            elseif subdir || in(name, dirs)
                # resumable functions doesn't support recursive generators
                for f in collect(walkfiles(path; exts, ex_dirs, dirs, subdir=true))
                    @yield f
                end
            end
            # file, only included
        elseif in(extension(splitpath(p)[end]), exts)
            @yield path
        end
    end
end


text:
"""
    GridData <: StaticArray

Simulation data specific to a single grid.

These behave like arrays, but contain both source and
destination arrays as simulations need separate read and
write steps to maintain independence between cells.

`GridData` objects also contain other data and settings needed
for optimisations.

# Type parameters

- `S`: grid size type tuple
- `R`: grid padding radius
- `T`: grid data type
"""
abstract type GridData{S,R,T,N} <: StaticArray{S,T,N} end

function (::Type{G})(d::GridData{S,R,T,N}) where {G<:GridData,S,R,T,N}
    args = source(d), dest(d), mask(d), proc(d), opt(d), boundary(d), padval(d), optdata(d)
    G{S,R,T,N,map(typeof, args)...}(args...)
end
function ConstructionBase.constructorof(::Type{T}) where T<:GridData{S,R} where {S,R}
    T.name.wrapper{S,R}
end

# Getters
radius(d::GridData{<:Any,R}) where R = R
mask(d::GridData) = d.mask
proc(d::GridData) = d.proc
opt(d::GridData) = d.opt
optdata(d::GridData) = d.optdata
boundary(d::GridData) = d.boundary
padval(d::GridData) = d.padval
source(d::GridData) = d.source
dest(d::GridData) = d.dest

# Get the size of the grid
gridsize(d::GridData) = size(d)
gridsize(A::AbstractArray) = size(A)
gridsize(nt::NamedTuple) = gridsize(first(nt))
gridsize(nt::NamedTuple{(),Tuple{}}) = 0, 0

# Get a view of the grid, without padding
gridview(d::GridData) = sourceview(d)
# Get a view of the grid source, without padding
sourceview(d::GridData) = _padless_view(source(d), axes(d), radius(d))
# Get a view of the grid dest, without padding
destview(d::GridData) = _padless_view(dest(d), axes(d), radius(d))

_padless_view(A::OffsetArray, axes, radius) = _padless_view(parent(A), axes, radius)
function _padless_view(A::AbstractArray, axes, radius)
    ranges = map(axes) do axis
        axis .+ radius
    end
    return view(A, ranges...)
end

# Get an a view of the source, preferring the underlying array if it is not a padded OffsetArray
source_array_or_view(d::GridData) = source(d) isa OffsetArray ? sourceview(d) : source(d)
# Get an a view of the dest, preferring the underlying array if it is not a padded OffsetArray
dest_array_or_view(d::GridData) = dest(d) isa OffsetArray ? destview(d) : dest(d)

# Base methods
function Base.copy!(grid::GridData{<:Any,R}, A::AbstractArray) where R
    pad_axes = map(ax -> ax .+ R, axes(A))
    copyto!(parent(source(grid)), CartesianIndices(pad_axes), A, CartesianIndices(A))
    return _update_optdata!(grid)
end
function Base.copy!(A::AbstractArray, grid::GridData{<:Any,R}) where R
    pad_axes = map(ax -> ax .+ R, axes(A))
    copyto!(A, CartesianIndices(A), parent(source(grid)), CartesianIndices(pad_axes))
    return A
end
function Base.copy!(A::AbstractDimArray{T,N}, grid::GridData{<:Any,R}) where {T,N,R}
    copy!(parent(A), grid)
    return A
end
function Base.copy!(grid::GridData{<:Any,R}, A::AbstractDimArray{T,N}) where {R,T,N}
    copy!(grid, parent(A))
    return grid
end
function Base.copy!(dst::GridData{<:Any,RD}, src::GridData{<:Any,RS}) where {RD,RS}
    dst_axes = map(s -> RD:s + RD, size(dst))
    src_axes = map(s -> RS:s + RS, size(src))
    copyto!(parent(source(dst)), CartesianIndices(dst_axes),
            parent(source(src)), CartesianIndices(src_axes)
    )
    return dst
end

@propagate_inbounds Base.getindex(d::GridData{s}, I...) where s = getindex(source(d), I...)
@propagate_inbounds function Base.getindex(d::GridData{s}, i1::Int, I::Int...) where s
    getindex(source(d), i1, I...)
end

# Local utility methods

# _addpadding => OffsetArray{T,N}
# Find the maximum radius required by all rules
# Add padding around the original init array, offset into the negative
# So that the first real cell is still 1, 1
function _addpadding(init::AbstractArray{T,1}, r, padval) where T
    l = length(init)
    paddedsize = l + 2r
    paddedaxis = -r + 1:l + r
    sourceparent = fill(convert(T, padval), paddedsize)
    source = OffsetArray(sourceparent, paddedaxis)
    # Copy the init array to the middle section of the source array
    for i in 1:l
        @inbounds source[i] = init[i]
    end
    return source
end
function _addpadding(init::AbstractArray{T,2}, r, padval) where T
    h, w = size(init)
    paddedsize = h + 4r, w + 2r
    paddedaxes = -r + 1:h + 3r, -r + 1:w + r
    pv = convert(eltype(init), padval)
    sourceparent = similar(init, typeof(pv), paddedsize...)
    sourceparent .= Ref(pv)
    # Copy the init array to the middle section of the source array
    _padless_view(sourceparent, axes(init), r) .= init
    source = OffsetArray(sourceparent, paddedaxes...)
    return source
end

# _swapsource => ReadableGridData
# Swap source and dest arrays of a grid
_swapsource(d::Tuple) = map(_swapsource, d)
function _swapsource(grid::GridData)
    src = grid.source
    dst = grid.dest
    @set! grid.dest = src
    @set! grid.source = dst
    _swapoptdata(opt(grid), grid)
end

_swapoptdata(opt::PerformanceOpt, grid::GridData) = grid

_build_optdata(opt::PerformanceOpt, init, r) = nothing

_update_optdata!(grid::GridData) = _update_optdata!(grid, opt(grid))
_update_optdata!(grid, opt) = grid

# _indtoblock
# Convert regular index to block index
@inline _indtoblock(x::Int, blocksize::Int) = (x - 1) ÷ blocksize + 1

# _blocktoind
# Convert block index to regular index
@inline _blocktoind(x::Int, blocksize::Int) = (x - 1) * blocksize + 1

"""
    ReadableGridData <: GridData

    ReadableGridData(grid::GridData)
    ReadableGridData{S,R}(init::AbstractArray, mask, opt, boundary, padval)

[`GridData`](@ref) object passed to rules for reading only.
Reads are always from the `source` array.
"""
struct ReadableGridData{
    S<:Tuple,R,T,N,Sc,D,M,P<:Processor,Op<:PerformanceOpt,Bo,PV,OD
} <: GridData{S,R,T,N}
    source::Sc
    dest::D
    mask::M
    proc::P
    opt::Op
    boundary::Bo
    padval::PV
    optdata::OD
end
function ReadableGridData{S,R}(
    source::Sc, dest::D, mask::M, proc::P, opt::Op, boundary::Bo,
    padval::PV, optdata::OD
) where {S,R,Sc<:AbstractArray{T,N},D,M,P,Op,Bo,PV,OD} where {T,N}
    ReadableGridData{S,R,T,N,Sc,D,M,P,Op,Bo,PV,OD}(
        source, dest, mask, proc, opt, boundary, padval, optdata
    )
end
@inline function ReadableGridData{S,R}(
    init::AbstractArray, mask, proc, opt, boundary, padval
) where {S,R}
    # If the grid radius is larger than zero we pad it as an OffsetArray
    if R > 0
        source = _addpadding(init, R, padval)
        dest = _addpadding(init, R, padval)
    else
        source = deepcopy(init)
        dest = deepcopy(init)
    end
    optdata = _build_optdata(opt, source, R)

    grid = ReadableGridData{S,R}(
        source, dest, mask, proc, opt, boundary, padval, optdata
    )
    return _update_optdata!(grid)
end

function Base.parent(d::ReadableGridData{S,<:Any,T,N}) where {S,T,N}
    SizedArray{S,T,N}(source_array_or_view(d))
end

"""
    WritableGridData <: GridData

    WritableGridData(grid::GridData)

[`GridData`](@ref) object passed to rules as write grids, and can be written
to directly as an array, or preferably using `add!` etc. All writes handle
updates to `SparseOpt()` and writing to the correct source/dest array.

Reads are always from the `source` array, while writes are always to the
`dest` array. This is because rules application must not be sequential
between cells - the order of cells the rule is applied to does not matter.
This means that using e.g. `+=` is not supported. Instead use `add!`.
"""
struct WritableGridData{
    S<:Tuple,R,T,N,Sc<:AbstractArray{T,N},D<:AbstractArray{T,N},
    M,P<:Processor,Op<:PerformanceOpt,Bo,PV,OD
} <: GridData{S,R,T,N}
    source::Sc
    dest::D
    mask::M
    proc::P
    opt::Op
    boundary::Bo
    padval::PV
    optdata::OD
end
function WritableGridData{S,R}(
    source::Sc, dest::D, mask::M, proc::P, opt::Op,
    boundary::Bo, padval::PV, optdata::OD
) where {S,R,Sc<:AbstractArray{T,N},D<:AbstractArray{T,N},M,P,Op,Bo,PV,OD} where {T,N}
    WritableGridData{S,R,T,N,Sc,D,M,P,Op,Bo,PV,OD}(
        source, dest, mask, proc, opt, boundary, padval, optdata
    )
end

Base.parent(d::WritableGridData{S}) where {S} = SizedArray{S}(destview(d))

### UNSAFE / LOCKS required

# Base.setindex!
# This is not safe for general use.
# It can be used where only identical transformations of a cell
# can happen from any other cell, such as setting all 1s to 2.
@propagate_inbounds function Base.setindex!(d::WritableGridData, x, I...)
    _setindex!(d, proc(d), x, I...)
end
@propagate_inbounds function Base.setindex!(d::WritableGridData, x, i1::Int, I::Int...)
    _setindex!(d, proc(d), x, i1, I...)
end

@propagate_inbounds function _setindex!(d::WritableGridData, proc::SingleCPU, x, I...)
    @boundscheck checkbounds(dest(d), I...)
    @inbounds _setoptindex!(d, x, I...)
    @inbounds dest(d)[I...] = x
end
@propagate_inbounds function _setindex!(d::WritableGridData, proc::ThreadedCPU, x, I...)
    # Dest status is not threadsafe, even if the
    # setindex itself is safe. So we LOCK
    lock(proc)
    @inbounds _setoptindex!(d, x, I...)
    unlock(proc)
    @inbounds dest(d)[I...] = x
end

# _setoptindex!
# Do anything the optimisation needs on `setindex`
_setoptindex!(d::WritableGridData{<:Any,R}, x, I...) where R = _setoptindex!(d, opt(d), x, I...)
_setoptindex!(d::WritableGridData{<:Any,R}, opt::PerformanceOpt, x, I...) where R = nothing


text:
plan_flows_fft(a::Array, args...; kwargs...) = plan_fft(a, args...; kwargs...)
plan_flows_rfft(a::Array, args...; kwargs...) = plan_rfft(a, args...; kwargs...)

"""
A placeholder grid object for `0D` problems (in other words, systems of ODEs).
"""
struct ZeroDGrid{T, A} <: AbstractGrid{T, A} end

"""
A one-dimensional `grid` object.
"""
struct OneDGrid{T<:AbstractFloat, Tk, Tx, Tfft, Trfft} <: AbstractGrid{T, Tk}
        nx :: Int
        nk :: Int
       nkr :: Int

        dx :: T
        Lx :: T

         x :: Tx

         k :: Tk
        kr :: Tk
    invksq :: Tk
   invkrsq :: Tk

   fftplan :: Tfft
  rfftplan :: Trfft

  # Range objects that access the aliased part of the wavenumber range
    kalias :: UnitRange{Int}
   kralias :: UnitRange{Int}
end

"""
    OneDGrid(nx, Lx; x0=-Lx/2, nthreads=Sys.CPU_THREADS, effort=FFTW.MEASURE,
                      T=Float64, dealias=1/3, ArrayType=Array)

Constructs a OneDGrid object with size `Lx`, resolution `nx`, and leftmost position `x0`.
FFT plans are generated for `nthreads` CPUs using FFTW flag `effort`. The float type is `T`
and the array types is `ArrayType`. The `dealias` keyword determines the highest  wavenubers
that are being zero-ed out by `dealias()` function; 1/3 is the nominal value for quadratic
nonlinearities.
"""
function OneDGrid(nx, Lx; x0=-Lx/2, nthreads=Sys.CPU_THREADS, effort=FFTW.MEASURE,
                  T=Float64, dealias=1/3, ArrayType=Array)

  dx = Lx/nx

   nk = nx
  nkr = Int(nx/2 + 1)

  # Physical grid
  x = range(T(x0), step=T(dx), length=nx)

  # Wavenubmer grid
   k = ArrayType{T}( fftfreq(nx, 2π/Lx*nx))
  kr = ArrayType{T}(rfftfreq(nx, 2π/Lx*nx))

   invksq = @. 1 / k^2
  invkrsq = @. 1 / kr^2
  CUDA.@allowscalar  invksq[1] = 0
  CUDA.@allowscalar invkrsq[1] = 0

  FFTW.set_num_threads(nthreads)
   fftplan = plan_flows_fft(ArrayType{Complex{T}, 1}(undef, nx), flags=effort)
  rfftplan = plan_flows_rfft(ArrayType{T, 1}(undef, nx), flags=effort)

  kalias, kralias = getaliasedwavenumbers(nk, nkr, dealias)

  Tx = typeof(x)
  Tk = typeof(k)
  Tfft = typeof(fftplan)
  Trfft = typeof(rfftplan)

  return OneDGrid{T, Tk, Tx, Tfft, Trfft}(nx, nk, nkr, dx, Lx, x, k, kr,
                                      invksq, invkrsq, fftplan, rfftplan, kalias, kralias)
end

"A two-dimensional `grid` object."
struct TwoDGrid{T<:AbstractFloat, Tk, Tx, Tfft, Trfft} <: AbstractGrid{T, Tk}
        nx :: Int
        ny :: Int
        nk :: Int
        nl :: Int
       nkr :: Int

        dx :: T
        dy :: T
        Lx :: T
        Ly :: T

         x :: Tx
         y :: Tx

         k :: Tk
         l :: Tk
        kr :: Tk
       Ksq :: Tk
    invKsq :: Tk
      Krsq :: Tk
   invKrsq :: Tk

   fftplan :: Tfft
  rfftplan :: Trfft

  # Range objects that access the aliased part of the wavenumber range
    kalias :: UnitRange{Int}
   kralias :: UnitRange{Int}
    lalias :: UnitRange{Int}
end

"""
    TwoDGrid(nx, Lx, ny=nx, Ly=Lx; x0=-Lx/2, y0=-Ly/2, nthreads=Sys.CPU_THREADS,
                      effort=FFTW.MEASURE, T=Float64, dealias=1/3, ArrayType=Array)

Constructs a TwoDGrid object with size `Lx`, `Ly`, resolution `nx`, `ny`, and leftmost
positions `x0`, `y0`. FFT plans are generated for `nthreads` CPUs using
FFTW flag `effort`. The float type is `T` and the array types is `ArrayType`.
The `dealias` keyword determines the highest wavenubers that are being zero-ed
out by `dealias()` function; 1/3 is the nominal value for quadratic nonlinearities.
"""
function TwoDGrid(nx, Lx, ny=nx, Ly=Lx; x0=-Lx/2, y0=-Ly/2, nthreads=Sys.CPU_THREADS,
                  effort=FFTW.MEASURE, T=Float64, dealias=1/3, ArrayType=Array)

  dx = Lx/nx
  dy = Ly/ny

   nk = nx
   nl = ny
  nkr = Int(nx/2 + 1)

  # Physical grid
  x = range(T(x0), step=T(dx), length=nx)
  y = range(T(y0), step=T(dy), length=ny)

  # Wavenubmer grid
   k = ArrayType{T}(reshape( fftfreq(nx, 2π/Lx*nx), (nk, 1)))
   l = ArrayType{T}(reshape( fftfreq(ny, 2π/Ly*ny), (1, nl)))
  kr = ArrayType{T}(reshape(rfftfreq(nx, 2π/Lx*nx), (nkr, 1)))

     Ksq = @. k^2 + l^2
  invKsq = @. 1 / Ksq
  CUDA.@allowscalar invKsq[1, 1] = 0

     Krsq = @. kr^2 + l^2
  invKrsq = @. 1 / Krsq
  CUDA.@allowscalar invKrsq[1, 1] = 0

  # FFT plans
  FFTW.set_num_threads(nthreads)
  fftplan = plan_flows_fft(ArrayType{Complex{T}, 2}(undef, nx, ny), flags=effort)
  rfftplan = plan_flows_rfft(ArrayType{T, 2}(undef, nx, ny), flags=effort)

  # Index endpoints for aliasfrac i, j wavenumbers
  kalias, kralias = getaliasedwavenumbers(nk, nkr, dealias)
  lalias, _ = getaliasedwavenumbers(nl, nl, dealias)

  Tx = typeof(x)
  Tk = typeof(k)
  Tfft = typeof(fftplan)
  Trfft = typeof(rfftplan)

  return TwoDGrid{T, Tk, Tx, Tfft, Trfft}(nx, ny, nk, nl, nkr, dx, dy, Lx, Ly, x, y, k, l, kr, Ksq, invKsq, Krsq, invKrsq,
                                      fftplan, rfftplan, kalias, kralias, lalias)
end

"A three-dimensional `grid` object."
struct ThreeDGrid{T<:AbstractFloat, Tk, Tx, Tfft, Trfft} <: AbstractGrid{T, Tk}
        nx :: Int
        ny :: Int
        nz :: Int
        nk :: Int
        nl :: Int
        nm :: Int
       nkr :: Int

        dx :: T
        dy :: T
        dz :: T
        Lx :: T
        Ly :: T
        Lz :: T

         x :: Tx
         y :: Tx
         z :: Tx

         k :: Tk
         l :: Tk
         m :: Tk
        kr :: Tk
       Ksq :: Tk
    invKsq :: Tk
      Krsq :: Tk
   invKrsq :: Tk

   fftplan :: Tfft
  rfftplan :: Trfft

  # Range objects that access the aliased part of the wavenumber range
    kalias :: UnitRange{Int}
   kralias :: UnitRange{Int}
    lalias :: UnitRange{Int}
    malias :: UnitRange{Int}
end

"""
    ThreeDGrid(nx, Lx, ny=nx, Ly=Lx, nz=nx, Lz=Lx; x0=-Lx/2, y0=-Ly/2, z0=-Lz/2,
    nthreads=Sys.CPU_THREADS, effort=FFTW.MEASURE, T=Float64, dealias=1/3, ArrayType=Array)

 Constructs a TwoDGrid object with size `Lx`, `Ly`, `Lz`, resolution `nx`, `ny`,
 `nz` and leftmost positions `x0`, `y0`, `z0`. FFT plans are generated for `nthreads`
 CPUs using FFTW flag `effort`. The float type is `T` and the array types is `ArrayType`.
 The `dealias` keyword determines the highest wavenubers that are being zero-ed
 out by `dealias()` function; 1/3 is the nominal value for quadratic nonlinearities.
"""
function ThreeDGrid(nx, Lx, ny=nx, Ly=Lx, nz=nx, Lz=Lx; x0=-Lx/2, y0=-Ly/2, z0=-Lz/2,
                  nthreads=Sys.CPU_THREADS, effort=FFTW.MEASURE, T=Float64, dealias=1/3, ArrayType=Array)

  dx = Lx/nx
  dy = Ly/ny
  dz = Lz/nz

  nk = nx
  nl = ny
  nm = nz
  nkr = Int(nx/2 + 1)

  # Physical grid
  x = range(T(x0), step=T(dx), length=nx)
  y = range(T(y0), step=T(dy), length=ny)
  z = range(T(z0), step=T(dz), length=nz)

  # Wavenubmer grid
   k = ArrayType{T}(reshape( fftfreq(nx, 2π/Lx*nx), (nk, 1, 1)))
   l = ArrayType{T}(reshape( fftfreq(ny, 2π/Ly*ny), (1, nl, 1)))
   m = ArrayType{T}(reshape( fftfreq(nz, 2π/Lz*nz), (1, 1, nm)))
  kr = ArrayType{T}(reshape(rfftfreq(nx, 2π/Lx*nx), (nkr, 1, 1)))

     Ksq = @. k^2 + l^2 + m^2
  invKsq = @. 1 / Ksq
  CUDA.@allowscalar invKsq[1, 1, 1] = 0

     Krsq = @. kr^2 + l^2 + m^2
  invKrsq = @. 1 / Krsq
  CUDA.@allowscalar invKrsq[1, 1, 1] = 0

  # FFT plans
  FFTW.set_num_threads(nthreads)
  fftplan = plan_flows_fft(ArrayType{Complex{T}, 3}(undef, nx, ny, nz), flags=effort)
  rfftplan = plan_flows_rfft(ArrayType{T, 3}(undef, nx, ny, nz), flags=effort)

  # Index endpoints for aliasfrac i, j wavenumbers
  kalias, kralias = getaliasedwavenumbers(nk, nkr, dealias)
  lalias, malias = getaliasedwavenumbers(nl, nm, dealias)

  Tx = typeof(x)
  Tk = typeof(k)
  Tfft = typeof(fftplan)
  Trfft = typeof(rfftplan)

  return ThreeDGrid{T, Tk, Tx, Tfft, Trfft}(nx, ny, nz, nk, nl, nm, nkr, dx, dy, dz, Lx, Ly, Lz,
                                        x, y, z, k, l, m, kr, Ksq, invKsq, Krsq, invKrsq, fftplan, rfftplan,
                                        kalias, kralias, lalias, malias)
end

Base.eltype(grid::OneDGrid) = eltype(grid.x)
Base.eltype(grid::TwoDGrid) = eltype(grid.x)
Base.eltype(grid::ThreeDGrid) = eltype(grid.x)

OneDGrid(dev::CPU, args...; kwargs...) = OneDGrid(args...; ArrayType=Array, kwargs...)
TwoDGrid(dev::CPU, args...; kwargs...) = TwoDGrid(args...; ArrayType=Array, kwargs...)
ThreeDGrid(dev::CPU, args...; kwargs...) = ThreeDGrid(args...; ArrayType=Array, kwargs...)

"""
    gridpoints(grid::TwoDGrid)
    gridpoints(grid::ThreeDGrid)

Returns the collocation points of the `grid` in 2D or 3D arrays `X, Y` (and `Z`).
"""
function gridpoints(grid::TwoDGrid{T, A}) where {T, A}
  X = [ grid.x[i₁] for i₁=1:grid.nx, i₂=1:grid.ny ]
  Y = [ grid.y[i₂] for i₁=1:grid.nx, i₂=1:grid.ny ]

  return A(X), A(Y)
end

function gridpoints(grid::ThreeDGrid{T, A}) where {T, A}
  X = [ grid.x[i₁] for i₁=1:grid.nx, i₂=1:grid.ny, i₃=1:grid.nz ]
  Y = [ grid.y[i₂] for i₁=1:grid.nx, i₂=1:grid.ny, i₃=1:grid.nz ]
  Z = [ grid.z[i₃] for i₁=1:grid.nx, i₂=1:grid.ny, i₃=1:grid.nz ]

  return A(X), A(Y), A(Z)
end

"""
    getaliasedwavenumbers(nk, nkr, aliasfraction)

Returns the top `aliasfraction` highest wavenumbers.
"""
function getaliasedwavenumbers(nk, nkr, aliasfraction)
  # Index endpoints for aliased i, j wavenumbers
  # 1/3 aliasfraction => upper 1/6 of +/- wavenumbers (1/3 total) are set to 0 after performing fft.
  # 1/2 aliasfraction => upper 1/4 of +/- wavenumbers (1/2 total) are set to 0 after performing fft.
  L = (1 - aliasfraction)/2 # (1 - 1/3) / 2 + 1 = 1/3.
  R = (1 + aliasfraction)/2 # (1 + 1/3) / 2 - 1 = 2/3.
  iL = floor(Int, L * nk) + 1
  iR =  ceil(Int, R * nk)

  aliasfraction < 1 || error("`aliasfraction` must be less than 1") # aliasfraction=1 is not sensible.
   kalias = (aliasfraction > 0) ? (iL:iR) : Int(nx/2 + 1)
  kralias = (aliasfraction > 0) ? (iL:nkr) : nkr

  return kalias, kralias
end

"""
    dealias!(a, g, kalias)

Dealias array `a` on the grid `g` with aliased x-wavenumbers `kalias`.
"""
function dealias!(a, g::OneDGrid)
  kalias = size(a, 1) == g.nkr ? g.kralias : g.kalias
  dealias!(a, g, kalias)
  return nothing
end

function dealias!(a, g::OneDGrid, kalias)
  @views @. a[kalias, :] = 0

  return nothing
end

function dealias!(a, g::TwoDGrid)
  kalias = size(a, 1) == g.nkr ? g.kralias : g.kalias
  dealias!(a, g, kalias)

  return nothing
end

function dealias!(a, g::TwoDGrid, kalias)
  @views @. a[kalias, g.lalias, :] = 0

  return nothing
end

function dealias!(a, g::ThreeDGrid)
  kalias = size(a, 1) == g.nkr ? g.kralias : g.kalias
  dealias!(a, g, kalias)

  return nothing
end

function dealias!(a, g::ThreeDGrid, kalias)
  @views @. a[kalias, g.lalias, g.malias, :] = 0

  return nothing
end

"""
    makefilter(K; order=4, innerK=0.65, outerK=1)

Returns a filter acting on the non-dimensional wavenumber K that decays exponentially
for K>innerK, thus removing high-wavenumber content from a spectrum it is multiplied
with. The decay rate is determined by order and outerK determines the outer wavenumber
at which the filter is smaller than Float64 machine precision.
"""
function makefilter(K::Array; order=4, innerK=0.65, outerK=1)
  TK = typeof(K)
  K = Array(K)
  decay = 15*log(10) / (outerK-innerK)^order # decay rate for filtering function
  filt = @. exp( -decay*(K-innerK)^order )
  filt[real.(K) .< innerK] .= 1

  return TK(filt)
end

function makefilter(g::OneDGrid; realvars=true, kwargs...)
  K = realvars ? g.kr*g.dx/π : @.(abs(g.k*g.dx/π))

  return makefilter(K; kwargs...)
end

function makefilter(g::TwoDGrid; realvars=true, kwargs...)
  K = realvars ?
      @.(sqrt((g.kr*g.dx/π)^2 + (g.l*g.dy/π)^2)) : @.(sqrt((g.k*g.dx/π)^2 + (g.l*g.dy/π)^2))

  return makefilter(K; kwargs...)
end

function makefilter(g::ThreeDGrid; realvars=true, kwargs...)
  K = realvars ?
      @.(sqrt((g.kr*g.dx/π)^2 + (g.l*g.dy/π)^2 + (g.m*g.dz/π)^2)) : @.(sqrt((g.k*g.dx/π)^2 + (g.l*g.dy/π)^2 + (g.m*g.dz/π)^2))

  return makefilter(K; kwargs...)
end

makefilter(g, T, sz; kwargs...) = ones(T, sz) .* makefilter(g; realvars=sz[1]==g.nkr, kwargs...)
makefilter(eq; kwargs...) = makefilter(eq.grid, fltype(eq.T), eq.dims; kwargs...)

"""
    griddevice(grid)

Returns the device on which the `grid` lives on.
"""
griddevice(grid::AbstractGrid{T, A}) where {T, A} = A<:Array ? "CPU" : "GPU"

show(io::IO, g::OneDGrid{T}) where T =
     print(io, "OneDimensionalGrid\n",
               "  ├─────────── Device: ", griddevice(g), '\n',
               "  ├──────── FloatType: $T", '\n',
               "  ├────────── size Lx: ", g.Lx, '\n',
               "  ├──── resolution nx: ", g.nx, '\n',
               "  ├── grid spacing dx: ", g.dx, '\n',
               "  └─────────── domain: x ∈ [$(g.x[1]), $(g.x[end])]")

show(io::IO, g::TwoDGrid{T}) where T =
     print(io, "TwoDimensionalGrid\n",
               "  ├───────────────── Device: ", griddevice(g), '\n',
               "  ├────────────── FloatType: $T", '\n',
               "  ├────────── size (Lx, Ly): ", (g.Lx, g.Ly), '\n',
               "  ├──── resolution (nx, ny): ", (g.nx, g.ny), '\n',
               "  ├── grid spacing (dx, dy): ", (g.dx, g.dy), '\n',
               "  └───────────────── domain: x ∈ [$(g.x[1]), $(g.x[end])]", '\n',
               "                             y ∈ [$(g.y[1]), $(g.y[end])]")

show(io::IO, g::ThreeDGrid{T}) where T =
     print(io, "ThreeDimensionalGrid\n",
               "  ├───────────────────── Device: ", griddevice(g), '\n',
               "  ├────────────────── FloatType: $T", '\n',
               "  ├────────── size (Lx, Ly, Lz): ", (g.Lx, g.Ly, g.Ly), '\n',
               "  ├──── resolution (nx, ny, nz): ", (g.nx, g.ny, g.nz), '\n',
               "  ├── grid spacing (dx, dy, dz): ", (g.dx, g.dy, g.dz), '\n',
               "  └────────────────────  domain: x ∈ [$(g.x[1]), $(g.x[end])]", '\n',
               "                                 y ∈ [$(g.y[1]), $(g.y[end])]", '\n',
               "                                 z ∈ [$(g.z[1]), $(g.z[end])]")


text:
using Statistics
using NCDatasets

using ArgParse
using JSON
using Formatting

include("LinearRegression.jl")
include("CESMReader.jl")
include("constants.jl")
using .CESMReader

correlation = (x1, x2) -> x1' * x2 / (sum(x1.^2)*sum(x2.^2)).^0.5

function parse_commandline()

    s = ArgParseSettings()
    @add_arg_table s begin

        "--data-file-prefix"
            help = "Data filename prefix including folder and path until the timestamp. File extension `nc` is assumed."
            arg_type = String
            required = true

        "--data-file-timestamp-form"
            help = "Data filename timestamp form. Either `YEAR` or `YEAR_MONTH`."
            arg_type = String
            required = true

        "--domain-file"
            help = "Domain file."
            arg_type = String
            required = true

        "--output-file"
            help = "Output atm temperature file."
            arg_type = String
            required = true

        "--beg-year"
            help = "Year of begin."
            arg_type = Int64
            required = true

        "--end-year"
            help = "Year of end."
            arg_type = Int64
            required = true

    end

    return parse_args(ARGS, s)
end

parsed = parse_commandline()
print(json(parsed, 4))

output_file = parsed["output-file"]

let

    if parsed["data-file-timestamp-form"] == "YEAR"
        filename_format = format("{:s}{{:04d}}.nc", joinpath(parsed["data-file-prefix"]))
        form = :YEAR
    elseif parsed["data-file-timestamp-form"] == "YEAR_MONTH"
        filename_format = format("{:s}{{:04d}}-{{:02d}}.nc", joinpath(parsed["data-file-prefix"]))
        form = :YEAR_MONTH
    end

    fh = FileHandler(filename_format=filename_format, form=form)

    #global hi   = getData(fh, "hi",   (parsed["beg-year"], parsed["end-year"]), (:, :))
    global aice = getData(fh, "ICEFRAC", (parsed["beg-year"], parsed["end-year"]), (:, :))

    #aice ./= 100.0

    global (Nx, Ny, Nt) = size(aice)
    global nyears = Int64(Nt / 12)

    if mod(Nt, 12) != 0
        ErrorException("Time record is not multiple of 12") |> throw
    end

end

Dataset(parsed["domain-file"], "r") do ds
    global mask = replace(ds["mask"], missing=>NaN)
    global area = replace(ds["area"], missing=>NaN)
    global lat  = replace(ds["yc"], missing=>NaN)

    global mask_idx = (mask .== 0.0)
    global total_area = sum(area)

    global NH_mask_idx = ( mask .== 0.0 ) .& ( lat .>= 0.0 )
    global SH_mask_idx = ( mask .== 0.0 ) .& ( lat .<  0.0 )

end

area_factor = 4π * Re^2.0 / total_area

total_ice_volume = zeros(Float64, Nt)
total_ice_area   = zeros(Float64, Nt)
total_ice_extent_area = zeros(Float64, Nt)

total_ice_volume_NH      = zeros(Float64, Nt)
total_ice_area_NH        = zeros(Float64, Nt)
total_ice_extent_area_NH = zeros(Float64, Nt)

total_ice_volume_SH      = zeros(Float64, Nt)
total_ice_area_SH        = zeros(Float64, Nt)
total_ice_extent_area_SH = zeros(Float64, Nt)

x = collect(Float64, 1:Nt)

# total ice volume and area
for i = 1:Nt
    a = area .* view(aice, :, :, i)
    #v = a .* view(hi, :, :, i)

    total_ice_area[i]     = sum( a[mask_idx]  )
    total_ice_area_NH[i]  = sum( a[NH_mask_idx]  )
    total_ice_area_SH[i]  = sum( a[SH_mask_idx]  )
    #total_ice_volume[i] = sum( tmp2[mask_idx] )
end

total_ice_area    .*= area_factor
total_ice_area_NH .*= area_factor
total_ice_area_SH .*= area_factor

# total ice extent area
for i = 1:Nt
    tmp = view(aice, :, :, i)

    extent_idx = tmp .>= .15

    total_ice_extent_area    = sum(area[mask_idx .& extent_idx])
    total_ice_extent_area_NH = sum(area[NH_mask_idx .& extent_idx]  )
    total_ice_extent_area_SH = sum(area[SH_mask_idx .& extent_idx]  )
end

total_ice_extent_area    .*= area_factor
total_ice_extent_area_NH .*= area_factor
total_ice_extent_area_SH .*= area_factor

Dataset(output_file, "c") do ds

    defDim(ds, "months", 12)
    defDim(ds, "time", Inf)
    defDim(ds, "Nx", Nx)
    defDim(ds, "Ny", Ny)

    for (varname, vardata, vardim, attrib) in [
        ("total_ice_volume",           total_ice_volume, ("time",), Dict()),
        ("total_ice_area",             total_ice_area  , ("time",), Dict()),
        ("total_ice_extent_area",      total_ice_extent_area  , ("time",), Dict()),
        ("total_ice_volume_NH",        total_ice_volume_NH, ("time",), Dict()),
        ("total_ice_volume_SH",        total_ice_volume_SH, ("time",), Dict()),
        ("total_ice_area_NH",          total_ice_area_NH  , ("time",), Dict()),
        ("total_ice_area_SH",          total_ice_area_SH  , ("time",), Dict()),
        ("total_ice_extent_area_NH",   total_ice_extent_area_NH  , ("time",), Dict()),
        ("total_ice_extent_area_SH",   total_ice_extent_area_SH  , ("time",), Dict()),
    ]

        println("Doing var: ", varname)

        var = defVar(ds, varname, Float64, vardim)
        var.attrib["_FillValue"] = 1e20

        var = ds[varname]

        for (k, v) in attrib
            var.attrib[k] = v
        end

        rng = []
        for i in 1:length(vardim)-1
            push!(rng, Colon())
        end
        push!(rng, 1:size(vardata)[end])
        var[rng...] = vardata

    end

end

max_stars_repo_name: meteorologytoday/CESM-diagnostic

max_issues_repo_name: meteorologytoday/CESM-diagnostic

max_forks_repo_name: meteorologytoday/CESM-diagnostic


text:
export JULIADYNAMICS_COLORS, JULIADYNAMICS_CMAP, JULIADYNAMICS_CMAP_DIVERGING
export lighten_color, darken_color, randomcolor

const MARKER = Circle(Point2f0(0, 0), Float32(1)) # allows pixel size (zoom independent)
const DEFAULT_BG = RGBf0(1.0, 1.0, 1.0)
using Makie: px

COLORSCHEME = [
    "#6F4AC7",
    "#2DB9C5",
    "#1B1B1B",
    "#E82727",
    "#004E41",
]
JULIADYNAMICS_COLORS = to_color.(COLORSCHEME)

"""
    darken_color(c, f = 1.2)
Darken given color `c` by a factor `f`.
If `f` is less than 1, the color is lightened instead.
"""
function darken_color(c, f = 1.2)
    c = to_color(c)
    return RGBAf0(clamp.((c.r/f, c.g/f, c.b/f, c.alpha), 0, 1)...)
end

"""
    lighten_color(c, f = 1.2)
Lighten given color `c` by a factor `f`.
If `f` is less than 1, the color is darkened instead.
"""
function lighten_color(c, f = 1.2)
    c = to_color(c)
    return RGBAf0(clamp.((c.r*f, c.g*f, c.b*f, c.alpha), 0, 1)...)
end

# JULIADYNAMICS_CMAP = [
#     lighten_color(JULIADYNAMICS_COLORS[3], 1.1),
#     darken_color(JULIADYNAMICS_COLORS[1], 1.2),
#     lighten_color(JULIADYNAMICS_COLORS[1], 1.2),
#     lighten_color(JULIADYNAMICS_COLORS[2], 1.0),
#     lighten_color(JULIADYNAMICS_COLORS[6], 3.5),
# ]
JULIADYNAMICS_CMAP = Reverse(:dense)

# JULIADYNAMICS_CMAP_DIVERGING = [
#     lighten_color(JULIADYNAMICS_COLORS[1], 2.2),
#     lighten_color(JULIADYNAMICS_COLORS[1], 1.0),
#     darken_color(JULIADYNAMICS_COLORS[3], 2.0),
#     lighten_color(JULIADYNAMICS_COLORS[2], 1.0),
#     lighten_color(JULIADYNAMICS_COLORS[2], 2.0),
# ]
JULIADYNAMICS_CMAP_DIVERGING = :curl


text:
function get_gps_orbital_elements(obj::GPSPropagator)
    return jcall(obj, "getGPSOrbitalElements", GPSOrbitalElements, ())
end


text:
module LibFaceDetection

using Colors
using FixedPointNumbers
using GeometryBasics
using libfacedetection_jll

# Write your package code here.
include("cwrapper.jl")

export detect_faces

end


text:
mutable struct Cvor
    color::Char
    type::Char
    d::Int
    pred::Int
end
mutable struct labyrinth
   rows::Int64
   cols::Int64
   AdjMatrix::Array{Int, 2}
   V::Array{Cvor, 1}
end

function generateTableAdjMatrix(nbRows, nbColumns)
    adjMatrix = zeros(nbRows*nbColumns, nbRows*nbColumns);
    adjMatrixRowCounter = 0;

    for i = 1:nbRows
       for j = 1:nbColumns
           currentId = (i-1)*nbColumns + j;
           neighbourRow = [currentId-1 currentId currentId+1];
           neighbourMatrix = [neighbourRow .- nbColumns; neighbourRow; neighbourRow .+ nbColumns];

           startRowIndex = 1;
           startColumnIndex = 1;
           endRowIndex = 3;
           endColumnIndex = 3;

           if i == 1
               startRowIndex = 2;
           end
           if j == 1
               startColumnIndex = 2;
           end
           if i == nbRows
              endRowIndex = 2;
           end
           if j == nbColumns
              endColumnIndex = 2;
           end

	   neighbours = [];
           for m = startRowIndex:endRowIndex
              for n = startColumnIndex:endColumnIndex
                  if neighbourMatrix[m,n]!=currentId
                    push!(neighbours, neighbourMatrix[m,n]);
                  end
              end
           end

           adjRow = zeros(1, nbRows*nbColumns);
           for m = 1:length(neighbours)
                adjRow[neighbours[m]] = 1;
           end

           adjMatrixRowCounter  = adjMatrixRowCounter + 1;
           adjMatrix[adjMatrixRowCounter, :] = adjRow;
       end
    end
    return adjMatrix
end

function emptyLabyrinth(nrow, ncol)
   AdjMatrix = generateTableAdjMatrix(nrow, ncol)
   V = Array{Cvor}(undef, nrow*ncol)
   for i in 1:nrow*ncol
     V[i] = Cvor('W', '*', 0, 0)
   end
   return labyrinth(nrow, ncol, AdjMatrix, V)
end

function printLabyrinthType(G)
  for i in 1:G.rows
    for j in 1:G.cols
      print(" $(G.V[(i-1)*G.cols+j].type)")
    end
  println()
  end
end

function printLabyrinthColor(G)
  for i in 1:G.rows
    for j in 1:G.cols
      print(" $(G.V[(i-1)*G.cols+j].color)")
    end
  println()
  end
end

function defineNodesType!(G, nodes, type )
    for u = nodes
        G.V[u].type = type;
    end
end

function printLabyrinthPath(G, path)
  for i in 1:G.rows
    for j in 1:G.cols
       if ((i-1)*G.cols+j) in path
         print(" *")
       else
         print(" $(G.V[(i-1)*G.cols+j].type)")
       end
    end
  println()
  end
end


text:
include("PolygonBase.jl")
include("silhouette&pieces.jl")
plotly()

# P = tri_l
# H = tri_s

display(note)

max_stars_repo_name: toshiki-kosen/SolvingTangramPuzzle

max_issues_repo_name: toshiki-kosen/SolvingTangramPuzzle

max_forks_repo_name: toshiki-kosen/SolvingTangramPuzzle


text:
"""
    base(v::AbstractVertex)

Return the vertex wrapped in `v` (if any).
"""
function base(::AbstractVertex) end

"""
    OutputsVertex

Decorates an AbstractVertex with output edges.
"""
struct OutputsVertex{V<:AbstractVertex} <: AbstractVertex
    base::V
    outs::Vector{AbstractVertex} # Untyped because we might add other vertices to it
end
OutputsVertex(v::AbstractVertex) = OutputsVertex(v, AbstractVertex[])
init!(v::OutputsVertex, p::AbstractVertex) = foreach(in -> push!(outputs(in), p), inputs(v))

base(v::OutputsVertex) = v.base
(v::OutputsVertex)(x...) = base(v)(x...)

inputs(v::OutputsVertex) = inputs(base(v))
outputs(v::OutputsVertex) = v.outs

# Must not copy the outs field as it is typically added by higher vertices
Functors.functor(::Type{<:OutputsVertex}, v) = (base = v.base,), newbase -> OutputsVertex(newbase[1])

"""
    InputSizeVertex

Vertex with an (immutable) size.

Intended use is for wrapping an InputVertex in conjuntion with mutation
"""
struct InputSizeVertex{V<:AbstractVertex} <: AbstractVertex
    base::V
    size::Int

    function InputSizeVertex(b::V, size::Integer) where V <:OutputsVertex
        this = new{V}(b, Int(size))
        init!(b, this)
        return this
    end
end
InputSizeVertex(name, size::Integer) = InputSizeVertex(InputVertex(name), size)
InputSizeVertex(b::AbstractVertex, size::Integer) = InputSizeVertex(OutputsVertex(b), size)

base(v::InputSizeVertex)::AbstractVertex = v.base
(v::InputSizeVertex)(x...) = base(v)(x...)

inputs(v::InputSizeVertex) = inputs(base(v))
outputs(v::InputSizeVertex) = outputs(base(v))

@functor InputSizeVertex

function Base.show(io::IO, v::InputSizeVertex)
    print(io, "InputSizeVertex(")
    show(io, base(v))
    print(io, ", ", v.size, ')')
end

"""
    MutationTrait

Base type for traits relevant when mutating.
"""
abstract type MutationTrait end

"""
    MutationSizeTrait
Base type for mutation traits relevant to size
"""
abstract type MutationSizeTrait <: MutationTrait end

"""
    SizeTransparent
Base type for mutation traits which are transparent w.r.t size, i.e size changes propagate both forwards and backwards.

Tip: Use with [`FixedSizeTrait`](@ref) if the function has parameters which must be aligned with the input and output sizes.
"""
abstract type SizeTransparent <: MutationSizeTrait end
"""
    SizeStack
Transparent size trait type where inputs are stacked, i.e output size is the sum of all input sizes.
"""
struct SizeStack <: SizeTransparent end
"""
    SizeInvariant
Transparent size trait type where all input sizes must be equal to the output size, e.g. elementwise operations (including broadcasted).
"""
struct SizeInvariant <: SizeTransparent end
"""
    SizeAbsorb
Size trait type for which size changes are absorbed, i.e they do not propagate forward.

Note that size changes do propagate backward as changing the input size of a vertex requires that the output size of its input is also changed and vice versa.
"""
struct SizeAbsorb <: MutationSizeTrait end
"""
    Immutable

Trait for vertices which are immutable. Typically inputs and outputs as those are fixed to the surroundings (e.g a data set).
"""
struct Immutable <: MutationTrait end

"""
    trait(v)

Return the [`MutationTrait`](@ref) for a vertex `v`.
"""
trait(::AbstractVertex) = Immutable()

"""
    DecoratingTrait <: MutationTrait

Avbstract trait which wraps another trait. The wrapped trait of a [`DecoratingTrait`](@ref) `t` is accessible through [`base(t)`](@ref base(t::DecoratingTrait)).
"""
abstract type DecoratingTrait <: MutationTrait end

"""
    base(t::DecoratingTrait)

Return the trait wrapped by `t`.
"""
base(t::DecoratingTrait) = t.base # Lets just guess if we end up here :)

"""
    NamedTrait <: DecoratingTrait
    NamedTrait(name, base)

Trait which attaches `name` to a vertex. Calling [`name(v)`](@ref) on a vertex with this trait returns `name`.
"""
struct NamedTrait{S, T<:MutationTrait} <: DecoratingTrait
    name::S
    base::T
end
base(t::NamedTrait) = t.base

name(t::DecoratingTrait) = name(base(t))
name(t::NamedTrait) = t.name
name(::MutationTrait) = nothing

@functor NamedTrait

function Base.show(io::IO, t::NamedTrait)
    print(io, "NamedTrait(")
    show(io, t.name)
    print(io, ", ")
    show(io, t.base)
    print(io, ')')
end

"""
    FixedSizeTrait <: DecoratingTrait

Trait which indicates that a vertex is [`SizeTransparent`](@ref) while still having a fixed size.

This prevents NaiveNASlib from inferring the size from neighbouring vertices.

As an example, the function `x -> 2 .* x` accepts any size of `x`, while the function `x -> [1,2,3] .* x`
is [`SizeInvariant`](@ref) but has a fixed size of 3.

Note that `FixedSizeTrait` does not imply that the vertex can't change size.
"""
struct FixedSizeTrait{T<:SizeTransparent} <: DecoratingTrait
    base::T
end
base(t::FixedSizeTrait) = t.base

"""
    AfterΔSizeTrait <: DecoratingTrait
    AfterΔSizeTrait(strategy::S, base::T)

Calls `after_Δnin(strategy, v, Δs, ischanged)` and `after_Δnout(strategy, v, Δ, ischanged)` after a size change for the vertex `v` which this trait is attached to.
"""
struct AfterΔSizeTrait{S, T<:MutationTrait} <: DecoratingTrait
    strategy::S
    base::T
end
base(t::AfterΔSizeTrait) = t.base

@functor AfterΔSizeTrait

"""
    MutationVertex

Vertex which may be subject to mutation.

The member trait describes the nature of the vertex itself, for example if size changes
are absorbed (e.g changing an `nin x nout` matrix to an `nin - Δ x nout` matrix) or if they
propagate to neighbouring vertices (and if so, how).
"""
struct MutationVertex{V<:AbstractVertex, T<:MutationTrait} <: AbstractVertex
    base::V
    trait::T

    function MutationVertex(b::V, t::T) where {V <: OutputsVertex, T <: MutationTrait}
        this = new{V, T}(b, t)
        init!(b, this)
        return this
    end
end
MutationVertex(b::AbstractVertex, t::MutationTrait) = MutationVertex(OutputsVertex(b), t)

@functor MutationVertex

base(v::MutationVertex) = v.base
trait(v::MutationVertex) = v.trait
(v::MutationVertex)(x...) = base(v)(x...)

inputs(v::MutationVertex)  = inputs(base(v))
outputs(v::MutationVertex) = outputs(base(v))

# Stuff for displaying information about vertices

show_less(io::IO, v::InputSizeVertex) = show_less(io, base(v))
show_less(io::IO, v::OutputsVertex) = show_less(io, base(v))

show_less(io::IO, v::MutationVertex) = show_less(io, trait(v), v)
show_less(io::IO, t::DecoratingTrait, v::AbstractVertex) = show_less(io, base(t), v)
show_less(io::IO, t::NamedTrait, ::AbstractVertex) = print(io, t.name)
show_less(io::IO, ::MutationTrait, v::AbstractVertex) = show_less(io, base(v))

function Base.show(io::IO, v::OutputsVertex; close=')')
     show(io, base(v); close="")
     print(io, ", outputs=")
     show(io, outputs(v))
     print(io, close)
 end

 function Base.show(io::IO, v::MutationVertex; close=')')
    print(io, "MutationVertex(")
    show(io, base(v))
    print(io, ", ")
    show(io, trait(v))
    print(io, close)
 end

# Stuff for logging

name(v::InputSizeVertex) = name(base(v))
name(v::OutputsVertex) = name(base(v))
name(v::MutationVertex) = name(trait(v), v)
name(t::MutationTrait, ::V) where V = string(nameof(V),  "::", summary(t))
name(t::NamedTrait, v) = t.name
name(t::DecoratingTrait, v) = name(base(t), v)

nameorrepr(v) = repr(v)
nameorrepr(v::InputSizeVertex) = nameorrepr(base(v))
nameorrepr(v::InputVertex) = name(v)
nameorrepr(v::MutationVertex) = nameorrepr(trait(v), v)
nameorrepr(t::DecoratingTrait, v) = nameorrepr(base(t), v)
nameorrepr(t::NamedTrait, v) = t.name
nameorrepr(::MutationTrait, v) = repr(v)

issizemutable(v::AbstractVertex) = issizemutable(trait(v))
issizemutable(t::DecoratingTrait) = issizemutable(base(t))
issizemutable(::MutationSizeTrait) = true
issizemutable(::Immutable) = false


text:
export nodecorrectness

# name => idx map like in statsbase
# indexmap(nodes) = Dict(node => i for (i,node) in enumerate(nodes))

"""
    nodecorrectness(f::AbstractVector{Int},
                    nodes1::AbstractVector,nodes2::AbstractVector) -> nc

Calculates node correctness when given an alignment.

# Arguments
- `f` : Alignment between nodes1 and nodes2. `f[i]` describes the aligned
node pairs `nodes1[i]` and `nodes2[f[i]]`. Thus, `f` describes `length(f)`
aligned node pairs.
- `nodes1`,`nodes2` : Node sets that `f` desribes the alignment of.

# Output
- `nc` : Node correctness between 0 and 1.
"""
nodecorrectness(f::AbstractVector{Int},
                nodes1::AbstractVector,nodes2::AbstractVector) =
                   mean(i -> nodes1[i]==nodes2[f[i]], 1:length(f))


text:
using Revise
using Images
using StaticArrays
using ChangePrecision
using BenchmarkTools
using Setfield
using CUDA

using RayTracingInOneJulia

import RayTracingInOneJulia.scatter

UnionMaterial = Union{Lambertian{T}, Metal{T}, Dielectric{T}} where T
struct UnionMaterialWrapper{T}; _::UnionMaterial{T}; end
scatter(r_in::Ray{T}, rec::HitRecord{T, UnionMaterialWrapper{T}}) where {T} = scatter(r_in, @set rec.material = rec.material._)

function main(use_cuda=true)
    @changeprecision Float32 begin
        T = typeof(0.0)

        # Image

        aspect_ratio = 16 / 9
        image_width = 400
        image_height = floor(Int, image_width / aspect_ratio)
        samples_per_pixel = 100
        max_depth = 50

        if use_cuda
            ArrType = CuArray
        else
            ArrType = Array
        end

        image = ArrType{RGB{T}}(undef, image_height, image_width)

        # World

        material_ground = Lambertian(RGB(0.8, 0.8, 0.0))
        material_center = Dielectric(1.5)
        material_left   = Dielectric(1.5)
        material_right  = Metal(RGB(0.8, 0.6, 0.2), 1.0)

        world = ArrType([
            Sphere(Point3( 0.0, -100.5, -1.0), 100.0, UnionMaterialWrapper(material_ground)),
            Sphere(Point3( 0.0,    0.0, -1.0),   0.5, UnionMaterialWrapper(material_center)),
            Sphere(Point3(-1.0,    0.0, -1.0),   0.5, UnionMaterialWrapper(material_left)),
            Sphere(Point3( 1.0,    0.0, -1.0),   0.5, UnionMaterialWrapper(material_right)),
            ])

        # Camera

        cam = Camera{T}()

        # Render

        @time begin
            render!(image, world, cam, image_width, image_height, samples_per_pixel, max_depth)
            ArrType==CuArray && CUDA.synchronize()
        end
        # @btime begin
        #     render!($image, $world, $cam, $image_width, $image_height, $samples_per_pixel, $max_depth)
        #     $ArrType==CuArray && CUDA.synchronize()
        # end

        image = ArrType != Array ? Array(image) : image

        # Save

        save("image.png", image)

        image

    end # @changeprecision
end

# #

main(false)
main(true)

# #

main(false)
main(false)
main(false)
main(true)
main(true)
main(true)


text:
using PulsarSearch
using Test
using Distributions
using Documenter

doctest(PulsarSearch)

@testset "PulsarSearch.jl" begin
    @testset "Basic" begin
        flatarr = Array([0.5])
        @test z_n(flatarr, 2) == 0
        flatarr = Array(range(1, stop = 10) / 10)
        @test z_n(flatarr, 2) < 1e-10
        @test z_n(ones(10), 2) == 40
    end

    @testset "Binned" begin
        flatarr = zeros(10)
        @test z_n_binned(flatarr, 2) == 0
        flatarr = ones(10)
        @test z_n_binned(flatarr, 2) < 1e-10
        arr = zeros(10)
        arr[1] = 10
        @test z_n(arr, 2) == 40
    end

    @testset "Search" begin
        phases = rand(Normal(0.5, 0.1), 10000)
        pulse_no = rand(Uniform(0, 1000), 10000)
        pulse_no = floor.(pulse_no)
        f = 1.123
        times = sort((phases + pulse_no) / f)
        freqs, stats = z_n_search(times, 2, 1.0, 1.5, oversample = 4.0)
        freqs_bin, stats_bin = z_n_search_hist(times, 2, 1.0, 1.5, oversample = 4.0, nbin=256)
        maxind = argmax(stats)
        @test abs(freqs[maxind] - f) < 1e-3
        maxind = argmax(stats_bin)
        @test abs(freqs_bin[maxind] - f) < 1e-3
        @test isapprox(stats_bin[maxind], stats[maxind], rtol=0.1)
    end
end

@testset "Stat" begin
    @testset "single_from_multi $ntrial" for ntrial in
                                             [1, 10, 100, 1000, 10000, 100000]
        epsilon_1 = 0.00000001
        epsilon_n = p_multitrial_from_single_trial(epsilon_1, ntrial)
        epsilon_1_corr = p_single_trial_from_p_multitrial(epsilon_n, ntrial)

        @test isapprox(epsilon_1_corr, epsilon_1; rtol = 1e-2)
    end
    @testset "Zn Det Lev" begin
        @test isapprox(z2_n_detection_level(2), 13.276704135987625)
        epsilon_corr = p_single_trial_from_p_multitrial(0.01, 2)
        @test isapprox(
            z2_n_detection_level(4, 0.01, ntrial = 2),
            z2_n_detection_level(4, epsilon_corr),
        )

    end
    @testset "Zn Probability ntrial $ntrial" for ntrial in
                                                 [1, 10, 100, 1000, 100000]
        detlev = z2_n_detection_level(2, 0.1, ntrial = ntrial)
        @test isapprox(z2_n_probability(detlev, 2, ntrial = ntrial), 0.1)
        @test isapprox(z2_n_logprobability(detlev, 2, ntrial = ntrial), -2.3025850929940455)
    end

    @testset "Chisq logp" begin
        chi2 = 31
        d = Chisq(2)
        @test isapprox(chi2_logp(chi2, 2), logccdf(d, chi2), atol = 0.1)
        chi2 = Array([5, 32])
        @test all(isapprox.(chi2_logp.(chi2, 2), logccdf.(d, chi2), atol = 0.1))
    end
end


text:
abstract hypobolicSolver_oneD <: oneDoperator

type upwind{n} <: hypobolicSolver_oneD
  f :: Function
end

type LaxWendroff <: hypobolicSolver_oneD
  f :: Function
end

type DiscreteData_oneD
  value :: Array{Float64,1}
  dx :: Float64
end

function eval!(p :: hypobolicSolver_oneD, data :: DiscreteData_oneD, t :: Float64, dt :: Float64)
  error("No Function for abstract types")
end

function eval(p :: upwind{1}, data :: DiscreteData_oneD , t :: Float64, dt :: Float64)
  y1 = data.value
  dx = data.dx
  fy1 = p.f(y1)
  fy2 = cat(1,fy1[end], fy1[1:end-1] )
  y2 = y1 - (fy1 - fy2) * dt/dx
  return DiscreteData_oneD(y2,dx)
end

function eval(p :: upwind{2}, data :: DiscreteData_oneD , t :: Float64, dt :: Float64)
  y1 = data.value
  dx = data.dx
  fy1 = p.f(y1)
  fy2 = cat(1, fy1[end], fy1[1:end-1])
  fy3 = cat(1, fy1[end-1:end], fy1[1:end-2])
  y2 = y1 - (+fy3 -4fy2 + 3fy1) * dt / (2dx)
  return DiscreteData_oneD(y2,dx)
end

function eval(p :: LaxWendroff, data :: DiscreteData_oneD , t :: Float64, dt :: Float64)
  y1 = data.value
  dx = data.dx
  fy1 = p.f(y1)

  yright = (cat(1, fy1[2:end], fy1[1])+ fy1)/2
  yleft = (cat(1, fy1[end], fy1[1:end-1])+ fy1)/2
  yright = yright - (cat(1, fy1[2:end], fy1[1]) - fy1) * dt / (2dx)
  yleft = yleft  - (fy1 - cat(1, fy1[end], fy1[1:end-1])) * dt / (2dx)

  y2 = y1 - (p.f(yright) - p.f(yleft)) * dt/dx
  return DiscreteData_oneD(y2,dx)
end


text:
# Note that the two calls to using Stats in this file are
# strictly required: one pulls the names into DataFrames
# for easy access within the module, whereas the other call
# pushes those names into Main.
using Stats
using DataArrays

module DataFrames

##############################################################################
##
## Dependencies
##
##############################################################################

using Base.Intrinsics
using DataArrays
using GZip
using Stats
using SortingAlgorithms

##############################################################################
##
## Global constants
##
##############################################################################

const DEFAULT_COLUMN_TYPE = Float64

##############################################################################
##
## Extend methods in Base by default
##
##############################################################################

importall Base
importall Stats
import Base: Sort, Order
import Base.Sort: sort, sort!, Algorithm, defalg, issorted
import Base.Order: Ordering, By, Lt, Perm, Forward, lt, ord
import Base.AsyncStream

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

export # reconcile_groups,
       @DataFrame,
       @transform,
       AbstractDataFrame,
       AbstractIndex,
       array,
       based_on,
       by,
       cbind,
       clean_colnames!,
       coefnames,
       colffts,
       colmaxs,
       colmeans,
       colmedians,
       colmins,
       colnorms,
       colprods,
       colranges,
       colstds,
       colsums,
       colvars,
       colwise,
       combine,
       complete_cases,
       complete_cases!,
       cut,
       DataFrame,
       DataStream,
       describe,
       dict,
       drop_duplicates!,
       duplicated,
       EachCol,
       EachRow,
       findat,
       flipud!,
       flipud,
       Formula,
       get_groups,
       gl,
       GroupApplied,
       groupby,
       GroupedDataFrame,
       Index,
       index,
       IndexedVector,
       Indexer,
       interaction_design_matrix,
       is_group,
       letters,
       LETTERS,
       load_df,
       matrix,
       merge,
       model_response,
       ModelFrame,
       ModelMatrix,
       nafilter,
       naFilter,
       NamedArray,
       names!,
       nareplace,
       naReplace,
       nas,
       ncol,
       nrow,
       order,
       paste,
       percent_change,
       pool,
       pool!,
       PooledDataVecs, # The capitalization and/or name for this is a bit inconsistent (merge_pools, maybe?). Do we want to export?
       printtable,
       range,
       rbind,
       read_minibatch,
       readtable,
       reldiff,
       rename_group!,
       rep,
       rename!,
       rename,
       rowffts,
       rowmaxs,
       rowmeans,
       rowmedians,
       rowmins,
       rownorms,
       rowprods,
       rowranges,
       rowstds,
       rowsums,
       rowvars,
       save,
       set_group,
       set_groups,
       SimpleIndex,
       stack,
       SubDataFrame,
       subset,
       table,
       types,
       unique,
       unstack,
       vector,
       with,
       within!,
       within,
       without,
       writetable,
       xtab,
       xtabs,
       stack_df,
       StackedVector,
       RepeatedVector,
       EachRepeatedVector,
       melt,
       melt_df,
       pivot_table,
       RComplex,     # Vector{Complex128}
       RInteger,     # Vector{Int32} plus BitVector of NA indicators
       RLogical,     # BitVector of values and BitVector of NA indicators
       RNumeric,     # Vector{Float64}
       RList,        # Vector{Any}
       RString,      # Vector{ASCIIString}
       RSymbol,      # Symbol stored as an String b/c of embedded '.'

       class,                              # in the S3 sense of "class"
       inherits,
       read_rda,
       vecbind

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

include("utils.jl")
include("index.jl")
include("namedarray.jl")
include("dataframe.jl")
include("show.jl")
include("merge.jl")
include("grouping.jl")
include("reshape.jl")
include("formula.jl")
include("io.jl")
include("datastream.jl")
include("operators.jl")
include("indexing.jl")
include("extras.jl")
include("RDA.jl")
include("dataframe_blocks.jl")

##############################################################################
##
## Deprecations
##
##############################################################################

Base.@deprecate vecbind vcat
Base.@deprecate rbind vcat
Base.@deprecate cbind hcat
Base.@deprecate read_table readtable
Base.@deprecate print_table printtable
Base.@deprecate write_table writetable
Base.@deprecate merge(df1::AbstractDataFrame, df2::AbstractDataFrame; on::Any = nothing, kind::Symbol = :inner) Base.join
Base.@deprecate colnames(adf::AbstractDataFrame) Base.names
Base.@deprecate colnames!(adf::AbstractDataFrame, vals) names!
Base.@deprecate coltypes(adf::AbstractDataFrame) types

end # module DataFrames


text:
# Don't paste into REPL, use `include("prompt_times.jl")`

using REPL
using Printf
@eval REPL begin
     const LAST_TIME = Ref{Float64}(time())
     const PROMPT_CACHE  = Ref{Union{String, Nothing}}(nothing)
end

# Format the time printing here as you wish
function prettytime(t)
    if REPL.PROMPT_CACHE[] !== nothing
        return REPL.PROMPT_CACHE[]
    end
    if t < 1e2
        value, units = t, "ns"
    elseif t < 1e5
        value, units = t / 1e3, "μs"
    elseif t < 1e8
        value, units = t / 1e6, "ms"
    elseif t < 120e9
        value, units = t / 1e9, "s"
    elseif t < 1e9*60*90
	value, units = t / 60e9, "min"
    else
	value, units = t / 3600e9, "hr"
    end
    s = string(round(value, sigdigits=3), " ", units)
    REPL.PROMPT_CACHE[] = s
    return s
end

julia_repl = Base.active_repl.interface.modes[1]
old_f = julia_repl.on_done
julia_repl.on_done = (args...) -> (REPL.LAST_TIME[] = time(); REPL.PROMPT_CACHE[] = nothing; old_f(args...))
julia_repl.prompt = () -> string("julia [", prettytime((time() - REPL.LAST_TIME[])*1e9), "]> ")


text:
using Markdown
using InteractiveUtils
using ReinforcementLearning
using Plots
using Flux
using Statistics

begin
	const NX = 4
	const NY = 12
	const Start = CartesianIndex(1, 1)
	const Goal = CartesianIndex(1, 12)
	const UDZ = [
        CartesianIndex(-1, 0),  # up
        CartesianIndex(1, 0),   # down
        CartesianIndex(0, 0),  # zero
	]
	const LinearInds = LinearIndices((NX, NY))
end

function iscliff(p::CartesianIndex{2})
    x, y = Tuple(p)
    x == 4 && y > 1 && y < NY
end

heatmap((!iscliff).(CartesianIndices((NX, NY))); yflip = true)

begin
	Base.@kwdef mutable struct CliffWalkingEnv <: AbstractEnv
		position::CartesianIndex{2} = Start
	end
	function (env::CliffWalkingEnv)(a::Int)
		x, y = Tuple(env.position + UDZ[a])
		env.position = CartesianIndex(min(max(x, 1), NX), min(max(y, 1), NY))
	end
end

RLBase.state(env::CliffWalkingEnv) = LinearInds[env.position]

RLBase.state_space(env::CliffWalkingEnv) = Base.OneTo(length(LinearInds))

RLBase.action_space(env::CliffWalkingEnv) = Base.OneTo(length(UDZ))

RLBase.reward(env::CliffWalkingEnv) = env.position == Goal ? 0.0 : (iscliff(env.position) ? -100.0 : -1.0)

RLBase.is_terminated(env::CliffWalkingEnv) = env.position == Goal || iscliff(env.position)

RLBase.reset!(env::CliffWalkingEnv) = env.position = Start

world = CliffWalkingEnv()

begin
	NS = length(state_space(world))
	NA = length(action_space(world))
end

create_agent(α, method) = Agent(
	policy = QBasedPolicy(
		learner=TDLearner(
			approximator=TabularQApproximator(
				;n_state=NS,
				n_action=NA,
				opt=Descent(α),
			),
			method=method,
			γ=1.0,
			n=0
		),
		explorer=EpsilonGreedyExplorer(0.1)
	),
	trajectory=VectorSARTTrajectory()
)

function repeated_run(α, method, N, n_episode, is_mean=true)
	env = CliffWalkingEnv()
	rewards = []
	for _ in 1:N
		h = TotalRewardPerEpisode(;is_display_on_exit=false)
		run(
			create_agent(α, method),
			env,
			StopAfterEpisode(n_episode;is_show_progress=false),
			h
		)
		push!(rewards, is_mean ? mean(h.rewards) : h.rewards)
	end
	mean(rewards)
end

α = 0.001
N = 500
n_episode = 500

begin
	p = plot(legend=:bottomright, xlabel="Episodes", ylabel="Sum of rewards during episode")
	plot!(p, repeated_run(α, :SARS, N, n_episode, false), label="QLearning")
	plot!(p, repeated_run(α, :SARSA, N, n_episode, false), label="SARSA")
	p
end

# begin
# 	A = 0.1:0.05:0.95
# 	fig_9_4 = plot(;legend=:bottomright, xlabel="α", ylabel="Sum of rewards per episode")

# 	plot!(fig_9_4, A, [repeated_run(α, :SARS, 100, 100) for α in A], linestyle=:dash ,markershape=:rect, label="Interim Q")
# 	plot!(fig_9_4, A, [repeated_run(α, :SARSA, 100, 100) for α in A], linestyle=:dash, markershape=:dtriangle, label="Interim SARSA")
# 	plot!(fig_9_4, A, [repeated_run(α, :ExpectedSARSA, 100, 100) for α in A], linestyle=:dash, markershape=:cross, label="Interim ExpectedSARSA")

# 	plot!(fig_9_4, A, [repeated_run(α, :SARS, 10, 5000) for α in A], linestyle=:solid ,markershape=:rect, label="Asymptotic interim Q")
# 	plot!(fig_9_4, A, [repeated_run(α, :SARSA, 10, 5000) for α in A], linestyle=:solid, markershape=:dtriangle, label="Asymptotic SARSA")
# 	plot!(fig_9_4, A, [repeated_run(α, :ExpectedSARSA, 10, 5000) for α in A], linestyle=:solid, markershape=:cross, label="Asymptotic ExpectedSARSA")
# 	fig_9_4
# end

max_stars_repo_name: cb-g/machine-learning-in-finance

max_issues_repo_name: cb-g/machine-learning-in-finance

max_forks_repo_name: cb-g/machine-learning-in-finance


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

# domains
include("plotrecipes/domains/point_set.jl")
include("plotrecipes/domains/regular_grid.jl")
include("plotrecipes/domains/structured_grid.jl")
include("plotrecipes/domains/view.jl")

# data
include("plotrecipes/data.jl")

# problems
include("plotrecipes/problems/estimation.jl")
include("plotrecipes/problems/simulation.jl")
include("plotrecipes/problems/learning.jl")

# solutions
include("plotrecipes/solutions/estimation.jl")
include("plotrecipes/solutions/simulation.jl")

# partitions and weights
include("plotrecipes/partitions.jl")
include("plotrecipes/weighting.jl")

# distribution plots
include("plotrecipes/distplot1D.jl")
include("plotrecipes/distplot2D.jl")
include("plotrecipes/cornerplot.jl")


text:
using LabelledArrays

export linearmap1

"""
    eom_linearmap1(x, p, n) → SVector{2}

# References
Chen, Yonghong, et al. "Analyzing multiple nonlinear time series with extended
Granger causality." Physics Letters A 324.1 (2004): 26-35
"""
function eom_linearmap1(x, p, t)
    c = p[1]
    x, y = (x...,)
    t = t + 3
    dx = 3.4*x*(t - 1)*(1 - x^2*(t - 1))*exp(-x^2*(t - 1)) + 0.8*x*(t - 2) + rand(Normal(0, 0.05))
    dy = 3.4*y*(t - 1)*(1 - y^2*(t - 1))*exp(-y^2*(t - 1)) + 0.5*y*(t - 2) + c*x*(t - 2) + rand(Normal(0, 0.05))
    return SVector{2}(dx, dy)
end

"""
    linearmap1(u₀, c) → DiscreteDynamicalSystem

# References
Chen, Yonghong, et al. "Analyzing multiple nonlinear time series with extended Granger causality." Physics Letters A 324.1 (2004): 26-35
"""
function linearmap1(u₀, c)
    p = @LArray [c] (:c)
    DiscreteDynamicalSystem(eom_linearmap1, u₀, p)
end

"""
    linearmap1(;u₀ = [1, rand(2)], c = 0.5) → DiscreteDynamicalSystem

# References
Chen, Yonghong, et al. "Analyzing multiple nonlinear time series with extended Granger causality." Physics Letters A 324.1 (2004): 26-35
"""
linearmap1(;u₀ = rand(2), c = 0.5) = linearmap1(u₀, c)


text:
module NCTiles

using NCDatasets,NetCDF,MeshArrays
using Dates,Printf,Pkg
#using MITgcmTools

version()=Pkg.TOML.parsefile(joinpath(dirname(pathof(NCTiles)), "..", "Project.toml"))["version"]

include("write.jl")
include("read.jl")
include("tile_support.jl")
include("helper_functions.jl")

export NCvar, BinData, NCData, TileData, readbin, readncfile
export createfile, addDim, addVar, addData, addDimData
export write

end # module


text:
function thermalize(filepath::String, r::MolecularDynamicsTrial, N::Int64,
        n::Int64 )
    for i in 1:n
        println("Thermalizing, step $i of $n, with $N verlet steps per relaxation")
        verlet!(r,N)
        save(filepath,r)
        relax!(r)
        save(filepath,r)
    end
end


text:
using Cairo

"""
    struct CairoScreen{S} <: AbstractScreen
A "screen" type for CairoMakie, which encodes a surface
and a context which are used to draw a Scene.
"""
struct CairoScreen{S}
    surface::S
    context::Cairo.CairoContext
end

function CairoScreen(w, h; device_scaling_factor = 1, antialias = Cairo.ANTIALIAS_BEST)
    w, h = (w, h) .* device_scaling_factor
    surf = Cairo.CairoARGBSurface(w, h)
    # this sets a scaling factor on the lowest level that is "hidden" so its even
    # enabled when the drawing space is reset for strokes
    # that means it can be used to increase or decrease the image resolution
    ccall((:cairo_surface_set_device_scale, Cairo.libcairo), Cvoid, (Ptr{Nothing}, Cdouble, Cdouble),
        surf.ptr, device_scaling_factor, device_scaling_factor)

    ctx = Cairo.CairoContext(surf)
    Cairo.set_antialias(ctx, antialias)

    return CairoScreen(surf, ctx)
end

function rgbatuple(c::Colorant)
    rgba = RGBA(c)
    red(rgba), green(rgba), blue(rgba), alpha(rgba)
end

to_2d_scale(x::Number) = Vec2f0(x)
to_2d_scale(x::Vec) = Vec2f0(x)

function project_position(scene, point, model)
    # use transform func
    res = scene.camera.resolution[]
    p4d = to_ndim(Vec4f0, to_ndim(Vec3f0, point, 0f0), 1f0)
    clip = scene.camera.projectionview[] * model * p4d
    @inbounds begin
    # between -1 and 1
        p = (clip ./ clip[4])[Vec(1, 2)]
        # flip y to match cairo
        p_yflip = Vec2f0(p[1], -p[2])
        # normalize to between 0 and 1
        p_0_to_1 = (p_yflip .+ 1f0) / 2f0
    end
    # multiply with scene resolution for final position
    return p_0_to_1 .* res
end

project_scale(scene, s::Number, model = Mat4f0(I)) = project_scale(scene, Vec2f0(s), model)

function project_scale(scene, s, model = Mat4f0(I))
    p4d = to_ndim(Vec4f0, s, 0f0)
    p = @inbounds (scene.camera.projectionview[] * model * p4d)[Vec(1, 2)] ./ 2f0
    return p .* scene.camera.resolution[]
end

function draw_marker(ctx, m, pos, scale, strokecolor, strokewidth, marker_offset)

    marker_offset = marker_offset + scale ./ 2

    pos += Point2f0(marker_offset[1], -marker_offset[2])

    # Cairo.scale(ctx, scale...)
    Cairo.move_to(ctx, pos[1] + scale[1]/2, pos[2])
    Cairo.arc(ctx, pos[1], pos[2], scale[1]/2, 0, 2*pi)
    Cairo.fill_preserve(ctx)

    Cairo.set_source_rgba(ctx, rgbatuple(strokecolor)...)
    Cairo.set_line_width(ctx, Float64(strokewidth))
    Cairo.stroke(ctx)
end

function draw_atomic(screen::CairoScreen, p::Scatter)
    ctx = screen.context
    model = p.transformation.model[]
    isempty(p.positions) && return
    size_model = p.transform_marker ? model : Mat4f0(I)

    # if we give size in pixels, the size is always equal to that value
    is_pixelspace = p.markerspace == Pixel
    broadcast_foreach(p.positions, p.color, p.markersize, p.strokecolor,
                      p.strokewidth, p.marker, p.markeroffset) do point, col,
                          markersize, strokecolor, strokewidth, marker, mo

        scale = if is_pixelspace
            to_2d_scale(markersize)
        else
            # otherwise calculate a scaled size
            project_scale(p, markersize, size_model)
        end
        offset = if is_pixelspace
            to_2d_scale(mo)
        else
            project_scale(p, mo, size_model)
        end

        pos = project_position(p, point, model)

        isnan(pos) && return

        Cairo.set_source_rgba(ctx, rgbatuple(col)...)

        draw_marker(ctx, marker, pos, scale, strokecolor, strokewidth, offset)
    end
    nothing
end


text:
using Documenter, ReverseDiff

makedocs(;
    modules=[ReverseDiff],
    sitename="ReverseDiff.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", nothing) == "true",
        canonical="http://www.juliadiff.org/ReverseDiff.jl",
    ),
    pages=[
        "Home" => "index.md",
        "Limitation of ReverseDiff" => "limits.md",
        "API" => "api.md",
    ],
    strict=true,
    checkdocs=:exports,
)

deploydocs(;
    repo="github.com/JuliaDiff/ReverseDiff.jl", push_preview=true,
)

max_stars_repo_name: jrevels/ReverseDiffPrototype

max_issues_repo_name: jrevels/ReverseDiffPrototype

max_forks_repo_name: jrevels/ReverseDiffPrototype


text:
export det

import ..DTree

using mlpack._Internal.io

import mlpack_jll
const detLibrary = mlpack_jll.libmlpack_julia_det

# Call the C binding of the mlpack det binding.
function det_mlpackMain()
  success = ccall((:det, detLibrary), Bool, ())
  if !success
    # Throw an exception---false means there was a C++ exception.
    throw(ErrorException("mlpack binding error; see output"))
  end
end

" Internal module to hold utility functions. "
module det_internal
  import ..detLibrary

import ...DTree

# Get the value of a model pointer parameter of type DTree.
function IOGetParamDTree(paramName::String)::DTree
  DTree(ccall((:IO_GetParamDTreePtr, detLibrary), Ptr{Nothing}, (Cstring,), paramName))
end

# Set the value of a model pointer parameter of type DTree.
function IOSetParamDTree(paramName::String, model::DTree)
  ccall((:IO_SetParamDTreePtr, detLibrary), Nothing, (Cstring, Ptr{Nothing}), paramName, model.ptr)
end

# Serialize a model to the given stream.
function serializeDTree(stream::IO, model::DTree)
  buf_len = UInt[0]
  buf_ptr = ccall((:SerializeDTreePtr, detLibrary), Ptr{UInt8}, (Ptr{Nothing}, Ptr{UInt}), model.ptr, Base.pointer(buf_len))
  buf = Base.unsafe_wrap(Vector{UInt8}, buf_ptr, buf_len[1]; own=true)
  write(stream, buf)
end
# Deserialize a model from the given stream.
function deserializeDTree(stream::IO)::DTree
  buffer = read(stream)
  DTree(ccall((:DeserializeDTreePtr, detLibrary), Ptr{Nothing}, (Ptr{UInt8}, UInt), Base.pointer(buffer), length(buffer)))
end
end # module

"""
    det(; [folds, input_model, max_leaf_size, min_leaf_size, path_format, skip_pruning, test, training, verbose])

This program performs a number of functions related to Density Estimation Trees.
 The optimal Density Estimation Tree (DET) can be trained on a set of data
(specified by `training`) using cross-validation (with number of folds specified
with the `folds` parameter).  This trained density estimation tree may then be
saved with the `output_model` output parameter.

The variable importances (that is, the feature importance values for each
dimension) may be saved with the `vi` output parameter, and the density
estimates for each training point may be saved with the `training_set_estimates`
output parameter.

Enabling path printing for each node outputs the path from the root node to a
leaf for each entry in the test set, or training set (if a test set is not
provided).  Strings like 'LRLRLR' (indicating that traversal went to the left
child, then the right child, then the left child, and so forth) will be output.
If 'lr-id' or 'id-lr' are given as the `path_format` parameter, then the ID
(tag) of every node along the path will be printed after or before the L or R
character indicating the direction of traversal, respectively.

This program also can provide density estimates for a set of test points,
specified in the `test` parameter.  The density estimation tree used for this
task will be the tree that was trained on the given training points, or a tree
given as the parameter `input_model`.  The density estimates for the test points
may be saved using the `test_set_estimates` output parameter.

# Arguments

 - `folds::Int`: The number of folds of cross-validation to perform for
      the estimation (0 is LOOCV)  Default value `10`.

 - `input_model::DTree`: Trained density estimation tree to load.
 - `max_leaf_size::Int`: The maximum size of a leaf in the unpruned, fully
      grown DET.  Default value `10`.

 - `min_leaf_size::Int`: The minimum size of a leaf in the unpruned, fully
      grown DET.  Default value `5`.

 - `path_format::String`: The format of path printing: 'lr', 'id-lr', or
      'lr-id'.  Default value `lr`.

 - `skip_pruning::Bool`: Whether to bypass the pruning process and output
      the unpruned tree only.  Default value `false`.

 - `test::Array{Float64, 2}`: A set of test points to estimate the density
      of.
 - `training::Array{Float64, 2}`: The data set on which to build a density
      estimation tree.
 - `verbose::Bool`: Display informational messages and the full list of
      parameters and timers at the end of execution.  Default value `false`.

# Return values

 - `output_model::DTree`: Output to save trained density estimation tree
      to.
 - `tag_counters_file::String`: The file to output the number of points
      that went to each leaf.  Default value ``.

 - `tag_file::String`: The file to output the tags (and possibly paths)
      for each sample in the test set.  Default value ``.

 - `test_set_estimates::Array{Float64, 2}`: The output estimates on the
      test set from the final optimally pruned tree.
 - `training_set_estimates::Array{Float64, 2}`: The output density
      estimates on the training set from the final optimally pruned tree.
 - `vi::Array{Float64, 2}`: The output variable importance values for each
      feature.

"""
function det(;
             folds::Union{Int, Missing} = missing,
             input_model::Union{DTree, Missing} = missing,
             max_leaf_size::Union{Int, Missing} = missing,
             min_leaf_size::Union{Int, Missing} = missing,
             path_format::Union{String, Missing} = missing,
             skip_pruning::Union{Bool, Missing} = missing,
             test = missing,
             training = missing,
             verbose::Union{Bool, Missing} = missing,
             points_are_rows::Bool = true)
  # Force the symbols to load.
  ccall((:loadSymbols, detLibrary), Nothing, ());

  IORestoreSettings("Density Estimation With Density Estimation Trees")

  # Process each input argument before calling mlpackMain().
  if !ismissing(folds)
    IOSetParam("folds", convert(Int, folds))
  end
  if !ismissing(input_model)
    det_internal.IOSetParamDTree("input_model", convert(DTree, input_model))
  end
  if !ismissing(max_leaf_size)
    IOSetParam("max_leaf_size", convert(Int, max_leaf_size))
  end
  if !ismissing(min_leaf_size)
    IOSetParam("min_leaf_size", convert(Int, min_leaf_size))
  end
  if !ismissing(path_format)
    IOSetParam("path_format", convert(String, path_format))
  end
  if !ismissing(skip_pruning)
    IOSetParam("skip_pruning", convert(Bool, skip_pruning))
  end
  if !ismissing(test)
    IOSetParamMat("test", test, points_are_rows)
  end
  if !ismissing(training)
    IOSetParamMat("training", training, points_are_rows)
  end
  if verbose !== nothing && verbose === true
    IOEnableVerbose()
  else
    IODisableVerbose()
  end

  IOSetPassed("output_model")
  IOSetPassed("tag_counters_file")
  IOSetPassed("tag_file")
  IOSetPassed("test_set_estimates")
  IOSetPassed("training_set_estimates")
  IOSetPassed("vi")
  # Call the program.
  det_mlpackMain()

  return det_internal.IOGetParamDTree("output_model"),
         Base.unsafe_string(IOGetParamString("tag_counters_file")),
         Base.unsafe_string(IOGetParamString("tag_file")),
         IOGetParamMat("test_set_estimates", points_are_rows),
         IOGetParamMat("training_set_estimates", points_are_rows),
         IOGetParamMat("vi", points_are_rows)
end


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

mutable struct ApplicationGatewayWebApplicationFirewallConfiguration <: SwaggerModel
    enabled::Any # spec type: Union{ Nothing, Bool } # spec name: enabled
    firewallMode::Any # spec type: Union{ Nothing, String } # spec name: firewallMode
    ruleSetType::Any # spec type: Union{ Nothing, String } # spec name: ruleSetType
    ruleSetVersion::Any # spec type: Union{ Nothing, String } # spec name: ruleSetVersion
    disabledRuleGroups::Any # spec type: Union{ Nothing, Vector{ApplicationGatewayFirewallDisabledRuleGroup} } # spec name: disabledRuleGroups
    requestBodyCheck::Any # spec type: Union{ Nothing, Bool } # spec name: requestBodyCheck
    maxRequestBodySize::Any # spec type: Union{ Nothing, Int32 } # spec name: maxRequestBodySize
    maxRequestBodySizeInKb::Any # spec type: Union{ Nothing, Int32 } # spec name: maxRequestBodySizeInKb
    fileUploadLimitInMb::Any # spec type: Union{ Nothing, Int32 } # spec name: fileUploadLimitInMb
    exclusions::Any # spec type: Union{ Nothing, Vector{ApplicationGatewayFirewallExclusion} } # spec name: exclusions

    function ApplicationGatewayWebApplicationFirewallConfiguration(;enabled=nothing, firewallMode=nothing, ruleSetType=nothing, ruleSetVersion=nothing, disabledRuleGroups=nothing, requestBodyCheck=nothing, maxRequestBodySize=nothing, maxRequestBodySizeInKb=nothing, fileUploadLimitInMb=nothing, exclusions=nothing)
        o = new()
        validate_property(ApplicationGatewayWebApplicationFirewallConfiguration, Symbol("enabled"), enabled)
        setfield!(o, Symbol("enabled"), enabled)
        validate_property(ApplicationGatewayWebApplicationFirewallConfiguration, Symbol("firewallMode"), firewallMode)
        setfield!(o, Symbol("firewallMode"), firewallMode)
        validate_property(ApplicationGatewayWebApplicationFirewallConfiguration, Symbol("ruleSetType"), ruleSetType)
        setfield!(o, Symbol("ruleSetType"), ruleSetType)
        validate_property(ApplicationGatewayWebApplicationFirewallConfiguration, Symbol("ruleSetVersion"), ruleSetVersion)
        setfield!(o, Symbol("ruleSetVersion"), ruleSetVersion)
        validate_property(ApplicationGatewayWebApplicationFirewallConfiguration, Symbol("disabledRuleGroups"), disabledRuleGroups)
        setfield!(o, Symbol("disabledRuleGroups"), disabledRuleGroups)
        validate_property(ApplicationGatewayWebApplicationFirewallConfiguration, Symbol("requestBodyCheck"), requestBodyCheck)
        setfield!(o, Symbol("requestBodyCheck"), requestBodyCheck)
        validate_property(ApplicationGatewayWebApplicationFirewallConfiguration, Symbol("maxRequestBodySize"), maxRequestBodySize)
        setfield!(o, Symbol("maxRequestBodySize"), maxRequestBodySize)
        validate_property(ApplicationGatewayWebApplicationFirewallConfiguration, Symbol("maxRequestBodySizeInKb"), maxRequestBodySizeInKb)
        setfield!(o, Symbol("maxRequestBodySizeInKb"), maxRequestBodySizeInKb)
        validate_property(ApplicationGatewayWebApplicationFirewallConfiguration, Symbol("fileUploadLimitInMb"), fileUploadLimitInMb)
        setfield!(o, Symbol("fileUploadLimitInMb"), fileUploadLimitInMb)
        validate_property(ApplicationGatewayWebApplicationFirewallConfiguration, Symbol("exclusions"), exclusions)
        setfield!(o, Symbol("exclusions"), exclusions)
        o
    end
end # type ApplicationGatewayWebApplicationFirewallConfiguration

const _property_map_ApplicationGatewayWebApplicationFirewallConfiguration = Dict{Symbol,Symbol}(Symbol("enabled")=>Symbol("enabled"), Symbol("firewallMode")=>Symbol("firewallMode"), Symbol("ruleSetType")=>Symbol("ruleSetType"), Symbol("ruleSetVersion")=>Symbol("ruleSetVersion"), Symbol("disabledRuleGroups")=>Symbol("disabledRuleGroups"), Symbol("requestBodyCheck")=>Symbol("requestBodyCheck"), Symbol("maxRequestBodySize")=>Symbol("maxRequestBodySize"), Symbol("maxRequestBodySizeInKb")=>Symbol("maxRequestBodySizeInKb"), Symbol("fileUploadLimitInMb")=>Symbol("fileUploadLimitInMb"), Symbol("exclusions")=>Symbol("exclusions"))
const _property_types_ApplicationGatewayWebApplicationFirewallConfiguration = Dict{Symbol,String}(Symbol("enabled")=>"Bool", Symbol("firewallMode")=>"String", Symbol("ruleSetType")=>"String", Symbol("ruleSetVersion")=>"String", Symbol("disabledRuleGroups")=>"Vector{ApplicationGatewayFirewallDisabledRuleGroup}", Symbol("requestBodyCheck")=>"Bool", Symbol("maxRequestBodySize")=>"Int32", Symbol("maxRequestBodySizeInKb")=>"Int32", Symbol("fileUploadLimitInMb")=>"Int32", Symbol("exclusions")=>"Vector{ApplicationGatewayFirewallExclusion}")
Base.propertynames(::Type{ ApplicationGatewayWebApplicationFirewallConfiguration }) = collect(keys(_property_map_ApplicationGatewayWebApplicationFirewallConfiguration))
Swagger.property_type(::Type{ ApplicationGatewayWebApplicationFirewallConfiguration }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_ApplicationGatewayWebApplicationFirewallConfiguration[name]))}
Swagger.field_name(::Type{ ApplicationGatewayWebApplicationFirewallConfiguration }, property_name::Symbol) =  _property_map_ApplicationGatewayWebApplicationFirewallConfiguration[property_name]

const _allowed_ApplicationGatewayWebApplicationFirewallConfiguration_firewallMode = ["Detection", "Prevention"]

function check_required(o::ApplicationGatewayWebApplicationFirewallConfiguration)
    (getproperty(o, Symbol("enabled")) === nothing) && (return false)
    (getproperty(o, Symbol("ruleSetType")) === nothing) && (return false)
    (getproperty(o, Symbol("ruleSetVersion")) === nothing) && (return false)
    true
end

function validate_property(::Type{ ApplicationGatewayWebApplicationFirewallConfiguration }, name::Symbol, val)
    if name === Symbol("firewallMode")
        Swagger.validate_param(name, "ApplicationGatewayWebApplicationFirewallConfiguration", :enum, val, _allowed_ApplicationGatewayWebApplicationFirewallConfiguration_firewallMode)
    end
    if name === Symbol("maxRequestBodySize")
        Swagger.validate_param(name, "ApplicationGatewayWebApplicationFirewallConfiguration", :maximum, val, 128, false)
        Swagger.validate_param(name, "ApplicationGatewayWebApplicationFirewallConfiguration", :minimum, val, 8, false)
    end
    if name === Symbol("maxRequestBodySizeInKb")
        Swagger.validate_param(name, "ApplicationGatewayWebApplicationFirewallConfiguration", :maximum, val, 128, false)
        Swagger.validate_param(name, "ApplicationGatewayWebApplicationFirewallConfiguration", :minimum, val, 8, false)
    end
    if name === Symbol("fileUploadLimitInMb")
        Swagger.validate_param(name, "ApplicationGatewayWebApplicationFirewallConfiguration", :minimum, val, 0, false)
    end
end


text:
using Molecules
using Test

@testset "Molecules" begin
    include("test_parsing.jl")
    include("test_transformations.jl")
    include("test_properties.jl")
    include("test_sym.jl")
end


text:
using Printf  # for Julia version 1.0+

function isnarcissist(n, b=10)
    -1 < n || return false
    d = digits(n, base=b)
    m = length(d)
    n == mapreduce((x)->x^m, +, d)
end

function findnarcissist(verbose=false)
    goal = 25
    ncnt = 0
    verbose && println("Finding the first ", goal, " Narcissistic numbers:")
    for i in 0:typemax(1)
        isnarcissist(i) || continue
        ncnt += 1
        verbose && println(@sprintf "    %2d %7d" ncnt i)
        ncnt < goal || break
    end
end

findnarcissist()
@time findnarcissist(true)

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode