text:
# Helpers
function default_model_cache_location()
    return joinpath(pkgdir(DifferentiableStateSpaceModels), ".function_cache")
end

# TODO Add docstring
function make_perturbation_model(H; t, y, x, steady_states=nothing,
                                 steady_states_iv=nothing, Γ, Ω=nothing, η, Q=I, p=nothing,
                                 model_name,
                                 model_cache_location=default_model_cache_location(),
                                 overwrite_model_cache=false, verbose=true, max_order=2,
                                 save_ip=true, save_oop=false, # only does inplace by default
                                 skipzeros=true, fillzeros=false)
    @assert max_order ∈ [1, 2]
    @assert save_ip || save_oop

    module_cache_path = joinpath(model_cache_location, model_name * ".jl")

    # only load cache if the module isn't already loaded in memory
    if (isdefined(Main, Symbol(model_name)) && !overwrite_model_cache)
        verbose && printstyled("Using existing module $model_name\n"; color=:cyan)
        return module_cache_path = module_cache_path
    end

    # if path already exists
    if (ispath(module_cache_path) && !overwrite_model_cache)
        # path exists and not overwriting
        verbose &&
            printstyled("Model already generated at $module_cache_path\n"; color=:cyan)
        return module_cache_path
    end

    n_y = length(y)
    n_x = length(x)
    n = n_y + n_x
    n_p = length(p)
    @assert n_p > 0  # code written to have at least one parameter
    n_ϵ = size(η, 2)
    n_z = (Q == I) ? n : size(Q, 1)

    # TODO: error check that p, y, x has no overlap
    # Get the markovian variables and create substitutions
    y_subs = StructArray(make_substitutions.(t, y))
    x_subs = StructArray(make_substitutions.(t, x))
    y = y_subs.var
    x = x_subs.var
    y_p = y_subs.var_p
    x_p = x_subs.var_p
    y_ss = y_subs.var_ss
    x_ss = x_subs.var_ss
    subs = vcat(x_subs, y_subs)
    all_to_markov = vcat(subs.markov_t, subs.markov_tp1, subs.markov_inf)
    all_to_var = vcat(subs.tp1_to_var, subs.inf_to_var)
    # Helper to take the [z(∞) ~ expr] and become [z => expr] after substitutions
    function equations_to_dict(equations)
        return Dict(Symbol(substitute(substitute(eq.lhs, all_to_markov), all_to_var)) => Num(substitute(eq.rhs,
                                                                                all_to_markov))
                    for eq in equations)
    end

    # create functions in correct order
    ȳ = isnothing(steady_states) ? nothing :
        order_vector_by_symbols(equations_to_dict(steady_states), y_subs.symbol)
    x̄ = isnothing(steady_states) ? nothing :
         order_vector_by_symbols(equations_to_dict(steady_states), x_subs.symbol)
    ȳ_iv = isnothing(steady_states_iv) ? nothing :
           order_vector_by_symbols(equations_to_dict(steady_states_iv), y_subs.symbol)
    x̄_iv = isnothing(steady_states_iv) ? nothing :
            order_vector_by_symbols(equations_to_dict(steady_states_iv), x_subs.symbol)

    # Get any latex generated stuff we wish for pretty display of the model
    H_latex = latexify(H)
    steady_states_latex = latexify(steady_states)
    steady_states_iv_latex = latexify(steady_states_iv)

    # steady state requiers differentiation after substitution, and wrt [y; x]
    H = substitute.(H, Ref(all_to_markov))
    H̄ = deepcopy(H)
    H̄ = substitute_and_simplify(H̄, all_to_var)

    # Derivatives utilities return nothing if either argument nothing
    H̄_w = nested_differentiate(H̄, [y; x]) # differentiate post-substitution wrt w = [y;x], force a dense one
    H_yp = nested_differentiate(H, y_p)
    H_y = nested_differentiate(H, y)
    H_xp = nested_differentiate(H, x_p)
    H_x = nested_differentiate(H, x)
    Ψ = [Symbolics.hessian(f, [y_p; y; x_p; x]; simplify=true) for f in H]  # need for 1st order derivatives
    Ψ_yp = (max_order < 2) ? nothing : nested_differentiate(Ψ, y_p)
    Ψ_y = (max_order < 2) ? nothing : nested_differentiate(Ψ, y)
    Ψ_xp = (max_order < 2) ? nothing : nested_differentiate(Ψ, x_p)
    Ψ_x = (max_order < 2) ? nothing : nested_differentiate(Ψ, x)

    # The parameter derivatives are maps for dispatching by Symbol
    # utility function substitutes/simplifies because these aren't themselves differentiated
    function differentiate_to_dict(f, p)
        return Dict([Symbol(p_val) => substitute_and_simplify(nested_differentiate(f,
                                                                                p_val),
                                                              all_to_var) for p_val in p])
    end
    differentiate_to_dict(::Nothing, p) = nothing
    H_p = differentiate_to_dict(H, p)
    Γ_p = differentiate_to_dict(Γ, p)
    Ω_p = differentiate_to_dict(Ω, p)
    ȳ_p = differentiate_to_dict(ȳ, p)
    x̄_p = differentiate_to_dict(x̄, p)
    H_yp_p = differentiate_to_dict(H_yp, p)
    H_xp_p = differentiate_to_dict(H_xp, p)
    H_y_p = differentiate_to_dict(H_y, p)
    H_x_p = differentiate_to_dict(H_x, p)
    Ψ_p = (max_order < 2) ? nothing : differentiate_to_dict(Ψ, p)

    # apply substitutions and simplify.  nothing stays nothing
    H = substitute_and_simplify(H, all_to_markov)
    H_yp = substitute_and_simplify(H_yp, all_to_var)
    H_xp = substitute_and_simplify(H_xp, all_to_var)
    H_x = substitute_and_simplify(H_x, all_to_var)
    H_y = substitute_and_simplify(H_y, all_to_var)
    Ψ = substitute_and_simplify(Ψ, all_to_var)
    Ψ_yp = substitute_and_simplify(Ψ_yp, all_to_var)
    Ψ_y = substitute_and_simplify(Ψ_y, all_to_var)
    Ψ_xp = substitute_and_simplify(Ψ_xp, all_to_var)
    Ψ_x = substitute_and_simplify(Ψ_x, all_to_var)

    # Generate all functions, and rename using utility
    function build_named_function(f, name, args...; symbol_dispatch=nothing)
        expr = isnothing(symbol_dispatch) ?
               build_function(f, args...; skipzeros, fillzeros) :
               build_function(f, nothing, args...; skipzeros, fillzeros)
        ip_function = name_symbolics_function(expr[1], Symbol(name); inplace=false,
                                              symbol_dispatch)
        oop_function = name_symbolics_function(expr[2], Symbol(name * "!"); inplace=true,
                                               symbol_dispatch)
        return ip_function, oop_function
    end
    build_named_function(::Nothing, name, args...; symbol_dispatch=nothing) = nothing

    Γ_expr = build_named_function(Γ, "Γ", p)
    Ω_expr = build_named_function(Ω, "Ω", p)
    H_expr = build_named_function(H, "H", y_p, y, y_ss, x_p, x, x_ss, p)
    H_yp_expr = build_named_function(H_yp, "H_yp", y, x, p)
    H_y_expr = build_named_function(H_y, "H_y", y, x, p)
    H_xp_expr = build_named_function(H_xp, "H_xp", y, x, p)
    H_x_expr = build_named_function(H_x, "H_x", y, x, p)
    Ψ_expr = build_named_function(Ψ, "Ψ", y, x, p)
    Ψ_yp_expr = (max_order < 2) ? nothing : build_named_function(Ψ_yp, "Ψ_yp", y, x, p)
    Ψ_y_expr = (max_order < 2) ? nothing : build_named_function(Ψ_y, "Ψ_y", y, x, p)
    Ψ_xp_expr = (max_order < 2) ? nothing : build_named_function(Ψ_xp, "Ψ_xp", y, x, p)
    Ψ_x_expr = (max_order < 2) ? nothing : build_named_function(Ψ_x, "Ψ_x", y, x, p)
    H̄_expr = build_named_function(H̄, "H̄", [y; x], p)
    H̄_w_expr = build_named_function(H̄_w, "H̄_w", [y; x], p)
    ȳ_iv_expr = build_named_function(ȳ_iv, "ȳ_iv", p)
    x̄_iv_expr = build_named_function(x̄_iv, "x̄_iv", p)
    ȳ_expr = build_named_function(ȳ, "ȳ", p)
    x̄_expr = build_named_function(x̄, "x̄", p)

    # Derivatives dispatch by symbol, and contain dummy "nothing" argument to replace
    function build_function_to_dict(f, name, args...)
        return Dict(key => build_named_function(f[key], name, args...; symbol_dispatch=key)
                    for key in keys(f))
    end
    build_function_to_dict(::Nothing, name, args...) = nothing

    Γ_p_expr = build_function_to_dict(Γ_p, "Γ_p", p)
    Ω_p_expr = build_function_to_dict(Ω_p, "Ω_p", p)
    H_yp_p_expr = build_function_to_dict(H_yp_p, "H_yp_p", y, x, p)
    H_y_p_expr = build_function_to_dict(H_y_p, "H_y_p", y, x, p)
    H_xp_p_expr = build_function_to_dict(H_xp_p, "H_xp_p", y, x, p)
    H_x_p_expr = build_function_to_dict(H_x_p, "H_x_p", y, x, p)
    H_p_expr = build_function_to_dict(H_p, "H_p", y, x, p)
    ȳ_p_expr = build_function_to_dict(ȳ_p, "ȳ_p", p)
    x̄_p_expr = build_function_to_dict(x̄_p, "x̄_p", p)
    Ψ_p_expr = (max_order < 2) ? nothing : build_function_to_dict(Ψ_p, "Ψ_p", y, x, p)

    verbose && printstyled("Done Building Model\n"; color=:cyan)

    # Separate filenames for different orders and function types.  For example
    mkpath(model_cache_location)
    mkpath(joinpath(model_cache_location, model_name))
    # module_cache_path has the core module stuff and includes the others
    zero_order_oop_path = joinpath(model_cache_location, model_name, "zero_order_oop.jl")
    zero_order_ip_path = joinpath(model_cache_location, model_name, "zero_order_ip.jl")
    first_order_oop_path = joinpath(model_cache_location, model_name, "first_order_oop.jl")
    first_order_ip_path = joinpath(model_cache_location, model_name, "first_order_ip.jl")
    second_order_oop_path = joinpath(model_cache_location, model_name,
                                     "second_order_oop.jl")
    second_order_ip_path = joinpath(model_cache_location, model_name, "second_order_ip.jl")

    # Basic definitions are independent of the order
    open(module_cache_path, "w") do io
        write(io, "module $(model_name)\n")
        write(io, "using LinearAlgebra, SymbolicUtils, LaTeXStrings\n")  # SymbolicUtils used in the generated functions
        # TODO: Use Val{...} instead?  Useful for any dispatching?
        write(io, "const max_order = $max_order\n")
        write(io, "const n_y = $n_y\n")
        write(io, "const n_x = $n_x\n")
        write(io, "const n_p = $n_p\n")
        write(io, "const n_ϵ = $n_ϵ\n")
        write(io, "const n_z = $n_z\n")
        if n_ϵ == 1
            write(io, "const η = reshape($η, $n_x, $n_ϵ)\n")
        else
            write(io, "const η = $η\n")
        end
        write(io, "const Q = $Q\n")
        write(io, "const has_Ω = $(!isnothing(Ω))\n")
        write(io, "# Display definitions\n")
        write(io, "const x_symbols = $(x_subs.symbol)\n")
        write(io, "const y_symbols = $(y_subs.symbol)\n")
        write(io, "const u_symbols = $([y_subs.symbol; x_subs.symbol])\n")
        write(io, "const p_symbols = $(Symbol.(p))\n")
        write(io, "const H_latex = L\"$H_latex\"\n")
        write(io, "const steady_states_latex = L\"$steady_states_latex\"\n")
        write(io, "const steady_states_iv_latex = L\"$steady_states_iv_latex\"\n")

        write(io, "# Function definitions\n")
        save_oop && write(io, "include(\"$model_name/zero_order_oop.jl\")\n")
        save_ip && write(io, "include(\"$model_name/zero_order_ip.jl\")\n")
        save_oop && write(io, "include(\"$model_name/first_order_oop.jl\")\n")
        save_ip && write(io, "include(\"$model_name/first_order_ip.jl\")\n")
        save_oop &&
            max_order > 1 &&
            write(io, "include(\"$model_name/second_order_oop.jl\")\n")
        save_ip &&
            max_order > 1 &&
            write(io, "include(\"$model_name/second_order_ip.jl\")\n")
        return write(io, "end\n") # end module
    end

    # Zero order includes steady state calculations and derivatives
    save_ip && open(zero_order_ip_path, "w") do io
        write(io, string(Γ_expr[2]) * "\n\n")
        if isnothing(Ω)
            write(io, "const Ω! = nothing\n\n")
        else
            write(io, string(Ω_expr[2]) * "\n\n")
        end
        write(io, string(H̄_expr[2]) * "\n\n")
        write(io, string(H̄_w_expr[2]) * "\n\n")
        if isnothing(steady_states_iv)
            write(io, "ȳ_iv! = nothing\n\n")
            write(io, "x̄_iv! = nothing\n\n")
        else
            write(io, string(ȳ_iv_expr[2]) * "\n\n")
            write(io, string(x̄_iv_expr[2]) * "\n\n")
        end
        if isnothing(steady_states)
            write(io, "ȳ! = nothing\n\n")
            write(io, "x̄! = nothing\n\n")
            write(io, "ȳ_p! = nothing\n\n")
            write(io, "x̄_p! = nothing\n\n")
        else
            write(io, string(ȳ_expr[2]) * "\n\n")
            write(io, string(x̄_expr[2]) * "\n\n")
            foreach(fun -> write(io, string(fun[2]) * "\n\n"), values(ȳ_p_expr))
            foreach(fun -> write(io, string(fun[2]) * "\n\n"), values(x̄_p_expr))
        end
        write(io, "const steady_state! = nothing\n\n")
        foreach(fun -> write(io, string(fun[2]) * "\n\n"), values(Γ_p_expr))
        if isnothing(Ω)
            write(io, "const Ω_p = nothing\n\n")
        else
            foreach(fun -> write(io, string(fun[2]) * "\n\n"), values(Ω_p_expr))
        end
    end
    save_oop && open(zero_order_oop_path, "w") do io
        write(io, string(Γ_expr[1]) * "\n\n")
        if isnothing(Ω)
            write(io, "const Ω = nothing\n\n")
        else
            write(io, string(Ω_expr[1]) * "\n\n")
        end
        write(io, string(H̄_expr[1]) * "\n\n")
        write(io, string(H̄_w_expr[1]) * "\n\n")
        write(io, string(ȳ_iv_expr[1]) * "\n\n")
        write(io, string(x̄_iv_expr[1]) * "\n\n")
        write(io, string(ȳ_expr[1]) * "\n\n")
        write(io, string(x̄_expr[1]) * "\n\n")
        write(io, "const steady_state = nothing\n\n")
        foreach(fun -> write(io, string(fun[1]) * "\n\n"), values(Γ_p_expr))
        if isnothing(Ω)
            write(io, "const Ω_p = nothing\n\n")
        else
            foreach(fun -> write(io, string(fun[1]) * "\n\n"), values(Ω_p_expr))
        end
        foreach(fun -> write(io, string(fun[1]) * "\n\n"), values(ȳ_p_expr))
        return foreach(fun -> write(io, string(fun[1]) * "\n\n"), values(x̄_p_expr))
    end

    # First order perturbations + d/dp
    save_ip && open(first_order_ip_path, "w") do io
        write(io, string(H_expr[2]) * "\n\n")
        write(io, string(H_yp_expr[2]) * "\n\n")
        write(io, string(H_y_expr[2]) * "\n\n")
        write(io, string(H_xp_expr[2]) * "\n\n")
        write(io, string(H_x_expr[2]) * "\n\n")
        write(io, string(Ψ_expr[2]) * "\n\n")
        foreach(fun -> write(io, string(fun[2]) * "\n\n"), values(H_yp_p_expr))
        foreach(fun -> write(io, string(fun[2]) * "\n\n"), values(H_y_p_expr))
        foreach(fun -> write(io, string(fun[2]) * "\n\n"), values(H_xp_p_expr))
        foreach(fun -> write(io, string(fun[2]) * "\n\n"), values(H_x_p_expr))
        return foreach(fun -> write(io, string(fun[2]) * "\n\n"), values(H_p_expr))
    end
    save_oop && open(first_order_oop_path, "w") do io
        write(io, string(H_expr[1]) * "\n\n")
        write(io, string(H_yp_expr[1]) * "\n\n")
        write(io, string(H_y_expr[1]) * "\n\n")
        write(io, string(H_xp_expr[1]) * "\n\n")
        write(io, string(H_x_expr[1]) * "\n\n")
        write(io, string(Ψ_expr[1]) * "\n\n")
        foreach(fun -> write(io, string(fun[1]) * "\n\n"), values(H_yp_p_expr))
        foreach(fun -> write(io, string(fun[1]) * "\n\n"), values(H_y_p_expr))
        foreach(fun -> write(io, string(fun[1]) * "\n\n"), values(H_xp_p_expr))
        foreach(fun -> write(io, string(fun[1]) * "\n\n"), values(H_x_p_expr))
        return foreach(fun -> write(io, string(fun[1]) * "\n\n"), values(H_p_expr))
    end

    # Second order perturbations + d/dp
    max_order > 1 &&
        save_ip &&
        open(second_order_ip_path, "w") do io
            write(io, string(Ψ_yp_expr[2]) * "\n\n")
            write(io, string(Ψ_y_expr[2]) * "\n\n")
            write(io, string(Ψ_xp_expr[2]) * "\n\n")
            write(io, string(Ψ_x_expr[2]) * "\n\n")
            return foreach(fun -> write(io, string(fun[2]) * "\n\n"), values(Ψ_p_expr))
        end

    max_order > 1 &&
        save_oop &&
        open(second_order_oop_path, "w") do io
            write(io, string(Ψ_yp_expr[1]) * "\n\n")
            write(io, string(Ψ_y_expr[1]) * "\n\n")
            write(io, string(Ψ_xp_expr[1]) * "\n\n")
            write(io, string(Ψ_x_expr[1]) * "\n\n")
            return foreach(fun -> write(io, string(fun[1]) * "\n\n"), values(Ψ_p_expr))
        end
    verbose && printstyled("Saved $model_name to $module_cache_path\n"; color=:cyan)

    # if the module was included, gets the module name, otherwise returns nothing
    return module_cache_path
end

# Utility for making model if required, loading it if already generated, etc.
macro make_and_include_perturbation_model(
    model_name, H, nt
)
    quote
        local model_name = $model_name
        local model_filename = "$(model_name).jl"
        local model_cache_path = joinpath(
            DifferentiableStateSpaceModels.default_model_cache_location(),
            model_filename,
        )
        if !isfile(model_cache_path)
            make_perturbation_model(
                H;
                model_name,
                $nt...
            )
            include(model_cache_path)
        elseif !isdefined(Main, Symbol(model_name))
            include(model_cache_path)
        end
        local mod = getfield(Main, Symbol(model_name))
        DifferentiableStateSpaceModels.PerturbationModel(mod)
    end |> esc
end


text:
"""
Associated with various objects that contain constant particle properties
"""
abstract type Particle end

"""
ID cosmic ray by number of protons, neutrons, and electrons
"""
struct ParticleID <: Particle
    p::Int64
    n::Int64
    e::Int64
end # ParticleID

fundamentalRestMassDict = Dict(
    :e => upreferred(511.0u"keV/c^2"),
    :p => upreferred(938.0u"MeV/c^2")
    )

"""
    mass(r::ParticleID; m::Dict = fundamentalRestMass)

Simple particle mass.

TODO: Account for binding energy
"""
mass(r::ParticleID; m::Dict = fundamentalRestMassDict) = r.p*m[:p]+r.e*m[:e] + 1.008665*m[:p]*r.n

"""
    charge(r::ParticleID)

Returns particle charge
"""
charge(r::ParticleID) = (r.p-r.e)*Unitful.q

Base.:(==)(a::ParticleID, b::ParticleID) = (a.p == b.p && a.n == b.n && a.e == b.e) ? true : false
Base.length(a::ParticleID) = 1 #defined for broadcasting purposes
Base.iterate(a::ParticleID) = (a,nothing)
export ParticleID, mass, charge

particle_dict = Dict(
                     "e-" => ParticleID(0,0,1),
                     "e+" => ParticleID(0,0,-1),
                     "p+" => ParticleID(1,0,0),
                     "p-" => ParticleID(-1,0,0),
                     "He+2" => ParticleID(2,2,0)
                    )
particle_keys = Vector(collect(keys(particle_dict)))
export particle_dict, particle_keys

ext: jl


text:
using StanSample, Test

ProjDir = @__DIR__
cd(ProjDir) do

  shared_file = "shared_funcs.stan"

  bernoulli_model = "
    functions{

      #include $(shared_file)
      //#include shared_funcs.stan // a comment
      //#include shared_funcs.stan // a comment
      //#include /Users/rob/shared_funcs.stan // a comment

      void model_specific_function(){
          real x = 1.0;
          return;
      }

    }
    data {
      int<lower=1> N;
      int<lower=0,upper=1> y[N];
    }
    parameters {
      real<lower=0,upper=1> theta;
    }
    model {
      model_specific_function();
      theta ~ beta(my_function(),1);
      y ~ bernoulli(theta);
    }
  ";

  stanmodel = SampleModel("parse_and_interpolate", bernoulli_model)

  observeddata = Dict("N" => 10, "y" => [0, 1, 0, 1, 0, 0, 0, 0, 0, 1])

  rc = stan_sample(stanmodel, data=observeddata)

  if success(rc)
    # Convert to an MCMCChains.Chains object
    samples = read_samples(stanmodel)

    # Show cmdstan summary in DataFrame format
    df = read_summary(stanmodel)

    # Retrieve mean value of theta from the summary
    @test df[df.parameters .== :theta, :mean][1] ≈ 0.33 rtol=0.1

  end

end


text:
####################
#
# Choropleth Plot
#
####################

# Plot choropleth  - given array of shapes and a canvas
function choropleth(shapes::AbstractArray{T, 1}, canvas::Compose.Context,
                    fill_data, fill_color_map;
                    convertcoords=lonlat_to_webmercator, img_width=12cm,
                    line_width=0.05mm, line_color="black",
                    transform=identity) where {T<:Shapefile.GeoInterface.AbstractGeometry}

    # This method is only for polygons
    @assert typeof(shapes[1]) <: ShpPolygon

    # Colours
    n_colors = length(fill_color_map)
    minval = transform(minimum(fill_data))
    maxval = transform(maximum(fill_data))
    range = maxval - minval

    # Plot outlines to canvas
    i = 1
    for polygon in shapes
        val = transform(Float64(isnan(fill_data[i]) ? 0.0 : fill_data[i]))  # CHECK CHANGE FROM ISNA TO ISNAN !!!
        i += 1
        normval = (val - minval) / range
        color_idx = round(Int, normval * (n_colors - 1) + 1)
        if length(polygon.parts) == 1
            canvas = compose_single(polygon, canvas, convertcoords, line_width, line_color, fill_color_map[color_idx])
        else
            canvas = compose_multi(polygon, canvas, convertcoords, line_width, line_color, fill_color_map[color_idx])
        end
    end

    # Return canvas
    canvas
end

# Plot choropleth - given array of shapes and an MBR
function choropleth(shapes::AbstractArray{T, 1}, MBR::Shapefile.Rect,
                    fill_data, fill_color_map;
                    convertcoords=lonlat_to_webmercator, img_width=12cm,
                    options...) where {T<:Shapefile.GeoInterface.AbstractGeometry}

    # Create canvas
    canvas = create_canvas(MBR, convertcoords, img_width)

    # Plot
    choropleth(shapes, canvas, fill_data, fill_color_map,
               convertcoords=convertcoords, img_width=img_width; options...)
end

# Plot choropleth - given an array of shapes
function choropleth(shapes::AbstractArray{T, 1}, fill_data, fill_color_map;
                    options...) where {T<:Shapefile.GeoInterface.AbstractGeometry}

    # Create MBR
    MBR = minBR(shapes)

    # Plot
    choropleth(shapes, MBR, fill_data, fill_color_map; options...)
end


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

using Markdown
using InteractiveUtils

# ╔═╡ 0d7e9e2a-9a1e-11ec-3f58-5938aac102fc
begin
    using Pkg
    Pkg.activate(; temp=true)
	packages = [
		PackageSpec(; name="CairoMakie", version="0.7")
		PackageSpec(; name="ColorTypes", version="0.11")
	]
    Pkg.add(packages)
	Pkg.develop(; path=dirname(dirname(dirname(@__DIR__))))
end

# ╔═╡ 59e351d0-8df4-4016-898e-3c4c6f554d17
using CairoMakie, ColorTypes, JModels

# ╔═╡ a88ccce3-34d3-483f-b313-595bab69ac07
md"""
This notebook generates the JModels logo and favicon.
Credits go to Makie.jl for making it easy to create a nice PNG with a transparent background.
"""

# ╔═╡ 91ae0221-d73e-44c5-a672-d7e360937e2e
julia_color = Dict(
	"blue" => RGB(77/255, 101/255, 175/255),
	"red" => RGB(201/255, 60/255, 50/255),
	"green" => RGB(59/255, 150/255, 71/255),
	"purple" => RGB(145/255, 89/255, 162/255)
)

# ╔═╡ 79c539d4-8788-43e0-bdd9-1065159a5a1d
size = 480;

# ╔═╡ 0990b969-b28d-45d2-89b1-2a1a8eb57ac4
logo = let
	fig = Figure(; resolution=(size, size), backgroundcolor=:transparent)
	ax = Axis(fig[1, 1]; backgroundcolor=:transparent)

	cluster(x::Int, y::Int) = Point2f.([x - 1, x, x + 1], [y, y + 2, y])

	markersize = 50

	scatter!(ax, cluster(2, 2); color=julia_color["red"], markersize)
	scatter!(ax, cluster(5, 4); color=julia_color["green"], markersize)
	scatter!(ax, cluster(8, 6); color=julia_color["purple"], markersize)

	limits!(ax, 0, 10, 1, 9)

	lines!(ax, [0.2, 9.8], [1.2, 8.8]; color=julia_color["blue"], linewidth=8)

	hidespines!(ax)
	hidedecorations!(ax)

	fig
end

# ╔═╡ 97d8d81c-0cf1-4102-ab83-0becfad30453
favicon_path = joinpath(pkgdir(JModels), "docs", "src", "assets", "favicon.png");

# ╔═╡ 3c461c35-d8d4-4dce-8b1b-276a64152b02
# Favicons are usually around 48x48.
let
	mkpath(dirname(favicon_path))
	save(favicon_path, logo, px_per_unit=48/size)
	cd(dirname(favicon_path)) do
		# Documenter.jl wants an ico file.
		cp("favicon.png", "favicon.ico"; force=true)
	end
end;

# ╔═╡ d0c67ad2-ef24-4405-b01a-b3f11e623383
logo_path = let
	path = joinpath(pkgdir(JModels), "docs", "src", "assets", "logo.png")
	mkpath(dirname(path))
	path
end

# ╔═╡ cd996076-d3c9-47e3-874a-fc70fa38ad86
save(logo_path, logo; px_per_unit=240/size);

# ╔═╡ Cell order:
# ╠═a88ccce3-34d3-483f-b313-595bab69ac07
# ╠═0d7e9e2a-9a1e-11ec-3f58-5938aac102fc
# ╠═59e351d0-8df4-4016-898e-3c4c6f554d17
# ╠═91ae0221-d73e-44c5-a672-d7e360937e2e
# ╠═79c539d4-8788-43e0-bdd9-1065159a5a1d
# ╠═0990b969-b28d-45d2-89b1-2a1a8eb57ac4
# ╠═97d8d81c-0cf1-4102-ab83-0becfad30453
# ╠═3c461c35-d8d4-4dce-8b1b-276a64152b02
# ╠═d0c67ad2-ef24-4405-b01a-b3f11e623383
# ╠═cd996076-d3c9-47e3-874a-fc70fa38ad86


text:
module RiskPerfTests

using Test
using Dates
using Statistics
using RiskPerf

# load test data produced using R PerformanceAnalytics package
asset_returns = parse.(Float64, readlines("data/recon_asset_rets.csv")[2:end])
market_returns = parse.(Float64, readlines("data/recon_market_rets.csv")[2:end])
rf = parse.(Float64, readlines("data/recon_rf_rets.csv")[2:end])

@test std(rf) ≈ 1.03135983021606e-07

@test skewness(asset_returns) ≈ -0.0904045109144944
@test skewness(asset_returns; method=:moment) ≈ -0.0904045109144944
@test skewness(asset_returns; method=:fisher_pearson) ≈ -0.09054037810903758                # bug in R's PerformanceAnalytics, corrected here
@test skewness(asset_returns; method=:sample) ≈ -0.0906763586376902

@test kurtosis(asset_returns) ≈ 0.0619856901585289
@test kurtosis(asset_returns; method=:excess) ≈ 0.0619856901585289
@test kurtosis(asset_returns; method=:moment) ≈ 3.06198569015853
@test kurtosis(asset_returns; method=:cornish_fisher) ≈ -0.158076535757453

@test downside_deviation(asset_returns, 0.001; method=:full) ≈ 0.00102022195386797
@test downside_deviation(asset_returns, rf; method=:full) ≈ 0.000474556578712896            # bug in R's PerformanceAnalytics, corrected here
@test downside_deviation(asset_returns, 0.001; method=:partial) ≈ 0.00123267724486424
@test downside_deviation(asset_returns, rf; method=:partial) ≈ 0.00085648312494845          # bug in R's PerformanceAnalytics, corrected here

@test lower_partial_moment(asset_returns, 0.001, 2, :full) ≈ 1.04085283515417e-06
@test lower_partial_moment(asset_returns, rf, 2, :full) ≈ 2.25203946399689e-07
@test lower_partial_moment(asset_returns, 0.001, 2, :partial) ≈ 1.51949319000609e-06
@test lower_partial_moment(asset_returns, rf, 2, :partial) ≈ 7.33563343321463e-07

@test higher_partial_moment(asset_returns, 0.001, 2, :full) ≈ 2.15200955473486e-07
@test higher_partial_moment(asset_returns, rf, 2, :full) ≈ 1.0554005124107e-06
@test higher_partial_moment(asset_returns, 0.001, 2, :partial) ≈ 6.83177636423766e-07
@test higher_partial_moment(asset_returns, rf, 2, :partial) ≈ 1.5229444623531e-06

@test upside_potential_ratio(asset_returns, 0.001; method=:full) ≈ 0.199853514985918        # bug in R's PerformanceAnalytics, corrected here
@test upside_potential_ratio(asset_returns, rf; method=:full) ≈ 1.49034594998309            # bug in R's PerformanceAnalytics, corrected here
@test upside_potential_ratio(asset_returns, 0.001; method=:partial) ≈ 0.525105446510415     # bug in R's PerformanceAnalytics, corrected here
@test upside_potential_ratio(asset_returns, rf; method=:partial) ≈ 1.19157955999556         # bug in R's PerformanceAnalytics, corrected here

@test volatility(asset_returns) ≈ 0.0010120151595378
@test tracking_error(asset_returns, market_returns) ≈ 0.00111905105002574
@test information_ratio(asset_returns, market_returns) ≈ 0.34034622188359
@test omega_ratio(asset_returns, 0.001) ≈ 0.297008426326568

@test sharpe_ratio(asset_returns) ≈ 0.511256629034303
@test sharpe_ratio(asset_returns; risk_free=rf) ≈ 0.501375462696111
@test adjusted_sharpe_ratio(asset_returns) ≈ 0.5069731154546914  # slightly different to R's PerformanceAnalytics due to their annualization (0.506476646827294)

@test sortino_ratio(asset_returns) ≈ 1.10000435875069
@test sortino_ratio(asset_returns; MAR=rf) ≈ 1.0692077438794

α1, β1 = capm(asset_returns, market_returns)
@test α1 ≈ 0.000512345109212185
@test β1 ≈ 0.0370188032088717

α2, β2 = capm(asset_returns, market_returns; risk_free=rf)
@test α2 ≈ 0.000502715256287843
@test β2 ≈ 0.037019963055605

γ = jensen_alpha(asset_returns, market_returns; risk_free=rf)
@test γ ≈ α2

γ = modified_jensen(asset_returns, market_returns; risk_free=rf)
@test γ ≈ α2/β2

@test treynor_ratio(asset_returns, market_returns) ≈ 0.013976666292467
@test treynor_ratio(asset_returns, market_returns; risk_free=rf) ≈ 0.0137061068404259

@test value_at_risk(asset_returns, 0.05, :historical) ≈ -0.00115939244815142
@test value_at_risk(asset_returns, 0.05, :gaussian) ≈ -0.00114638483011466
@test value_at_risk(asset_returns, 0.05, :cornish_fisher) ≈ -0.00117095813580772

@test expected_shortfall(asset_returns, 0.05, :historical) ≈ -0.00164648648602694
@test expected_shortfall(asset_returns, 0.05, :gaussian) ≈ -0.00156905316259027
@test expected_shortfall(asset_returns, 0.05, :cornish_fisher) ≈ -0.0016287316713825

@test length(drawdowns(asset_returns)) == length(asset_returns)
@test mean(drawdowns(asset_returns)) ≈ -0.000418043468523338
@test std(drawdowns(asset_returns)) ≈ 0.000757130672676854
@test mean(drawdowns(asset_returns; geometric=true)) ≈ -0.000517982256217951
@test std(drawdowns(asset_returns; geometric=true)) ≈ 0.000916577011725993

pnl = [-50, 100, -20, 30, -120, -10, 200, 10, 18, -12, -20, -30, 0];
@test length(drawdowns_pnl(pnl)) == length(pnl)
@test mean(drawdowns_pnl(pnl)) ≈ -37.53846153846154
@test std(drawdowns_pnl(pnl)) ≈ 45.40303841634277

# https://breakingdownfinance.com/finance-topics/modern-portfolio-theory/marginal-contribution-to-risk-mctr/
w = [0.042, 0.25, 0.32, 0.22, 0.168]
covmat = [
    0.001369 0.001184 0.000161 0.001924 0.002886;
    0.001184 0.006400 0.000580 0.008840 0.006240;
    0.000161 0.000580 0.000841 0.001885 0.001131;
    0.001924 0.008840 0.001885 0.016900 0.003380;
    0.002886 0.006240 0.001131 0.003380 0.067600
]
ctb = risk_contribution(w, covmat)

@test all(isapprox.(ctb, [0.009944 0.217649 0.059174 0.284805 0.428427]'; rtol=0.0001))

end


text:
using WaterModelsAnalytics

import Memento

const _WM = WaterModelsAnalytics._WM
const _IM = WaterModelsAnalytics._WM._IM
const _MOI = WaterModelsAnalytics._WM._MOI

# Suppress warnings during testing.
Memento.setlevel!(Memento.getlogger(_IM), "error")
Memento.setlevel!(Memento.getlogger(_WM), "error")
WaterModelsAnalytics.logger_config!("error")

using Test

# Setup common test data paths (from dependencies).
wm_path = joinpath(dirname(pathof(_WM)), "..")

@testset "WaterModelsAnalytics" begin

    include("graph.jl")

    include("io.jl")

end


text:
using VCFTools
using MendelImpute
using Random

"""
    filter_and_mask(data::String, samples::Int)

Creates reference haplotypes and (unphased) target genotype files from `data`.

# Inputs
`data`: The full (phased) data simulated by msprime.
`samples`: Number of samples (genotypes) desired in target file. Remaining haplotypes will become the reference panel
"""
function filter_and_mask(data::String, samples::String)
    missingprop = 0.1
    n = nsamples(data)
    p = nrecords(data)
    samples = convert
    samples > n && error("requested samples exceed total number of genotypes in $data.")

    # output filenames (tgt_data1.vcf.gz, ref_data1.vcf.gz, and tgt_masked_data1.vcf.gz)
    tgt = "./tgt_" * data
    ref = "./ref_" * data
    tgt_mask = "./tgt_masked_" * data
    tgt_mask_unphase = "./tgt_masked_unphased_" * data

    # compute target and reference index
    tgt_index = falses(n)
    tgt_index[1:samples] .= true
    ref_index = .!tgt_index
    record_index = 1:p # save all records (SNPs)

    # generate masking matrix with `missingprop`% of trues (true = convert to missing)
    Random.seed!(2020)
    masks = falses(p, samples)
    for j in 1:samples, i in 1:p
        rand() < missingprop && (masks[i, j] = true)
    end

    # create outputs
    VCFTools.filter(data, record_index, tgt_index, des = tgt)
    VCFTools.filter(data, record_index, ref_index, des = ref)
    mask_gt(tgt, masks, des=tgt_mask)

    # finally, unphase the target data
    unphase(tgt_mask, outfile=tgt_mask_unphase)
end

data = ARGS[1]
samples = parse(Int, ARGS[2])
filter_and_mask(data, samples)

max_stars_repo_name: Hua-Zhou/MendelImpute


text:
"""
    fitbayes(model::EVA; niter::Int=5000, warmup::Int=2000)::BayesianEVA

Fit the extreme value model under the Bayesian paradigm.

"""
function fitbayes(model::EVA; niter::Int=5000, warmup::Int=2000)::BayesianEVA

    # Set initial values to the maximum likelihood estimates
    ml = fit(model)
    initialvalues = ml.θ̂

    # Define the loglikelihood function and the gradient for the NUTS algorithm
    logf(θ::DenseVector) = loglike(model,θ)
    Δlogf(θ::DenseVector) = ForwardDiff.gradient(logf, θ)
    function logfgrad(θ::DenseVector)
        ll = logf(θ)
        g = Δlogf(θ)
        return ll, g
    end

    sim = Chains(niter, nparameter(model), start = (warmup + 1))
    θ = NUTSVariate(initialvalues, logfgrad)
    @showprogress for i in 1:niter
        sample!(θ, adapt = (i <= warmup))
        if i > warmup
            sim[i, :, 1] = θ
        end
    end

    fittedmodel = BayesianEVA(model, sim)

    return fittedmodel

end

include(joinpath("bayesian", "bayesian_gev.jl"))
include(joinpath("bayesian", "bayesian_gp.jl"))


text:
using Sundials, Plots
pyplot(size=(700,700))

ProjDir = dirname(@__FILE__)
cd(ProjDir) do

  ## f2 routine. Compute function f(x,y).

  function f2(x, y, ydot)
      ydot[1] = y[2]
      ydot[2] = 2.0*y[1] - 3.0*y[2] + 3.0*x^2
  end

  x = 0.0:0.2:5
  res = Sundials.cvode(f2, [1.0, 0.0], collect(x))
  res |> display
  println()

  p1 = plot(x, res[:,1], color=:darkblue,
    xlabel="x", ylabel="y[1]",
    title="Example Sundials01A (using Sundials.jl)"
    )
  p2 = plot(x, res[:,2], color=:darkred,
    xlabel="x", ylabel="y[2]"
    )
  plot(p1, p2, layout=(2, 1))
  savefig("Sundials01A.png")
  gui()

end


text:
module annotations

import SynapseClient: @annotationsfunction, @standalonefunction, @createtype, @dicttype, AbstractSynapseDict, synapsecall, synapseclient, pyannotations

using PyCall

export
	to_synapse_annotations,
	from_synapse_annotations,
	to_submission_status_annotations,
	from_submission_status_annotations,
	set_privacy

@dicttype Annotations pyannotations.Annotations

@annotationsfunction to_synapse_annotations
@annotationsfunction from_synapse_annotations
@annotationsfunction to_submission_status_annotations
@annotationsfunction from_submission_status_annotations
@annotationsfunction set_privacy

end


text:
module AlignSpectra

using HDF5, PyPlot, DynamicTimeWarp
using Dierckx, Optim

include("MonotoneSpline.jl")

export MonotoneSpline, MonotoneSplineLogLog, compose_spline,
    compose_splinelog, matchspectra

"""Find <npeaks> distinct peaks in a histogram with contents <h>.
The algorithm is to label the fullest bin as the 1st peak, and penalize
all other bins by P=d/(d+<distscale>), where d is the bin distance to the
nearest peak. The bin with the 2nd-highest value of P*h is the 2nd peak,
the penalties P are recomputed, and the selection step is repeated until
<npeaks> peaks are found and return as a sorted list of bin numbers.

Returns: peak_bin_numbers::Vector{eltype(h)}, a vector of length <npeaks>.
"""
function distinctpeaks(h::Vector, npeaks::Integer, distscale::Number)
    if npeaks > length(h)
        error("Histogram has only $(length(h)) bins: cannot find $(npeaks).")
    end
    bin = collect(1:length(h))
    peaks = Int[indmax(h)]
    dist = abs(peaks[1] .- bin)
    penalty = dist ./ (distscale+dist)
    for pknum in 2:npeaks
        for p in peaks[2:end]
            dist = abs(p .- bin)
            penalty = min(penalty, dist./(distscale+dist))
        end
        push!(peaks, indmax(h .* penalty))
    end
    sort(peaks)
end

"""Replace a DTW path (p1,p2) with (x,y) where the x values are unique and
the y values are an average over all those p2 where p1==x."""
function uniquepath(p1::Vector, p2::Vector)
    @assert length(p1)==length(p2)

    minx = p1[1]
    x = unique(p1)
    sort!(x)
    y = zeros(Float64, length(x))
    n = zeros(Int, length(x))

    bin = 1
    for i=1:length(p2)
        while x[bin] < p1[i]
            bin += 1
        end
        n[bin] += 1
        y[bin] += p2[i]
    end
    float(x), y./n
end

"""Find the best cubic spline with knots at {k}, such that the Euclidean distance to all pairs
of points {x_i, y_i} is minimized and such that the spline has "natural boundary conditions"
(i.e., it has zero second derivative at the extreme knots, k[1] and k[end]).
"""
function findbestspline(x::Vector, y::Vector, k::Vector)
    k = sort(k)
    const m = length(k)
    const n = length(x)
    @assert length(x) == length(y)

    Mdiag = ones(Float64, m)
    for i=2:m-1
        Mdiag[i] = 2.*(k[i+1]-k[i-1])
    end
    Mdu = zeros(Float64, m-1)
    Mdl = zeros(Float64, m-1)
    for i=1:m-2
        Mdu[i+1] = k[i+2]-k[i+1]
        Mdl[i] = k[i+1]-k[i]
    end
    M = Tridiagonal(Mdl, Mdiag, Mdu)

    Ndiag = zeros(Float64, m)
    for i=2:m-1
        Ndiag[i] = 6.*(1./(k[i]-k[i-1]) - 1./(k[i+1]-k[i]))
    end
    Ndl = zeros(Float64, m-1)
    Ndu = zeros(Float64, m-1)
    for i=2:m-1
        Ndl[i-1] = -6./(k[i]-k[i-1])
        Ndu[i] = +6./(k[i]-k[i-1])
    end
    N = Tridiagonal(Ndl, Ndiag, Ndu)

    # s'' = R s where s are the f(x) values at x={k} and s'' are 2nd derivatives there
    R = inv(M) * full(N)

    V = zeros(Float64, n, m)
    for i=1:n
        xi = x[i]

        # Below lowest knot
        if xi<=k[1]
            thisk = 1
            dr = k[2]-k[1]
            B=(xi-k[thisk])/dr
            A=1-B
            V[i,thisk] = A
            V[i,thisk+1] = B

        # Above highest knot
        elseif xi>=k[m]
            thisk = m-1
            dr = k[end]-k[end-1]
            B=(xi-k[thisk])/dr
            A=1-B
            V[i,thisk] = A
            V[i,thisk+1] = B

        # In the proper spline range
        else
            thisk = findn(xi.>=k)[1][1]
            if thisk>=m || thisk<1
                @show thisk, k, xi
            end
            dr = k[thisk+1]-k[thisk]
            B=(xi-k[thisk])/dr
            A=1-B
            C = (A^3-A)*(dr^2)/6.
            D = (B^3-B)*(dr^2)/6.
            V[i,thisk] = A
            V[i,thisk+1] = B
            V[i,:] += C*R[thisk,:]
            V[i,:] += D*R[thisk+1,:]
        end
    end
    knotvals = pinv(V)*y
end

"""Find the best cubic spline according to the Dierckx.Spline1D function with
weight vector w."""
function findbestspline_numeric(x::Vector, y::Vector, w::Vector)
    Dierckx.Spline1D(x, y, w=w, bc="nearest", s=1.0)
end

"""Find the best cubic spline with knots at {k}, according to Dierckx.Spline1D with
fixed knots.
"""
function findbestspline_knots(x::Vector, y::Vector, knots::Vector)
    @show knots
    Dierckx.Spline1D(x, y, knots, bc="nearest", k=3)
end

function basicdtw(data1::Vector{Float32}, data2::Vector{Float32}) #c1::Vector{Int}, c2::Vector{Int})
    clf()
    subplot(321)
    c1,b,_=plt.hist(data1, 8000, [0,8000], histtype="step", color="r", normed=true)
    db = 0.5*(b[2]-b[1])
    binctrs = b[1:end-1]+db
    pk1 = distinctpeaks(c1, 8, 100.)
    for p in pk1
        plt.plot(binctrs[p], c1[p], "ro")
    end

    c2,_,_=plt.hist(data2, 8000, [0,8000], histtype="step", color="k", normed=true)

    d1d2_distance = DynamicTimeWarp.Distance.poissonpenalty(length(data1), length(data2))
    dist,p1,p2 = DynamicTimeWarp.fastdtw(c1, c2, 10, d1d2_distance)
    subplot(322)
    u1,u2 = uniquepath(binctrs[p1],binctrs[p2])
    plot(u1,u2,"g")
    for p in pk1
        plt.plot(binctrs[p], u2[p], "ro")
    end

    f = splinelogspace(binctrs[pk1], float(u2[pk1]), bc="extrapolate")
    #   f = MonotoneSpline(binctrs[pk1], float(u2[pk1]), bc="extrapolate")
    plot(u1, f(u1), "k")

    data1e = f(float(data1))
    subplot(323)
    c1,b,_=plt.hist(data1e, 8000, [0,8000], histtype="step", color="r", normed=true)
    c2,_,_=plt.hist(data2, 8000, [0,8000], histtype="step", color="k", normed=true)
    dist,p1,p2 = DynamicTimeWarp.fastdtw(c1, c2, 20, d1d2_distance)
    subplot(324)
    u1,u2 = uniquepath(binctrs[p1],binctrs[p2])
    plot(u1,u2,"g")
    pk1 = distinctpeaks(c1, 8, 100.)
    for p in pk1
        plt.plot(binctrs[p], u2[p], "ro")
    end

    #   @show pk1, binctrs[pk1], u2[pk1]
    g = splinelogspace(binctrs[pk1], float(u2[pk1]), bc="extrapolate")
    plot(u1, g(u1), "k")

    data1f = g(float(data1e))
    subplot(313)
    c1,b,_=plt.hist(data1f, 8000, [0,8000], histtype="step", color="r", normed=true)
    c2,_,_=plt.hist(data2, 8000, [0,8000], histtype="step", color="k", normed=true)
    return data1f
end

"""Take a single DTW step to align data1 and data2.
Finds "distinct peaks" in data1 histogram (of a minimum size and separation).
Finds the FastDTW path to align the spectra (using a Poisson-inspired distance
function) via fastdtw and then "uniquifies" the path using uniquepath so that
no x value (data1 bin number) is hit twice. Finally, it returns only the
peak locations in the two spaces, but doing lookups in the unique path.

Return: (p1, p2) where p1 is a list of peak locations in the data1
histogram, and p2 is the location of these same peaks in data2, as
determined by FastDTW.

Uses a fixed histogram binning of 0:2:8000 at this point."""
function one_dtw_step(hist1::Vector, hist2::Vector, binctrs::Vector)

    pk1 = distinctpeaks(hist1, 8, 100.)

    d1d2_distance = DynamicTimeWarp.Distance.poissonpenalty(sum(hist1), sum(hist2))
    dist,p1,p2 = DynamicTimeWarp.fastdtw(hist1, hist2, 20, d1d2_distance)
    @printf("The Poisson DTW path distance: %f\n", dist)
    u1,u2 = uniquepath(binctrs[p1],binctrs[p2])
    return binctrs[pk1], float(u2[pk1])
end

function matchspectra{T<:Number}(histograms::Matrix{T}; histrange::Range=0:0.5:8000)
    N = size(histograms)[2]
    if N == 1
        dummyspline = AlignSpectra.MonotoneSplineLogLog([1,10], [1,10]; bc="extrapolate")
        return [dummyspline], histograms[:,1]
    end

    if N == 2
        binedges = collect(histrange)
        db = binedges[2]-binedges[1]
        binctrs = binedges[1:end-1]+db*0.5
        Nbins = length(binctrs)

        n1, n2 = sum(histograms[:,1]), sum(histograms[:,2])
        r,s = one_dtw_step(histograms[:,1], histograms[:,2], binctrs)
        middle = (r*n1 .+ s*n2) / (n1+n2)
        f1 = AlignSpectra.MonotoneSplineLogLog(r, middle; bc="extrapolate")
        f2 = AlignSpectra.MonotoneSplineLogLog(s, middle; bc="extrapolate")

        # Now convert histograms to equivalent values, transform them, and
        # re-histogram the result. Approximate this by random sampling for
        # allowed values in the bin (for bins with <10 samples), or by evenly
        # spreading out the values to
        combined = zeros(histograms[:,1])
        for i = 1:length(combined)
            nh = histograms[i,1]
            if nh > 0
                Nval = div(nh, 8) + zeros(Int, 8)
                for j=1:nh-sum(Nval)
                    Nval[rand(1:8)] += 1
                end
                @assert sum(Nval) == nh
                values = linspace(binedges[i]+db/16., binedges[i+1]-db/16, 8)
                fvalues = f1(values)
                for (n, fv) in zip(Nval, fvalues)
                    bnum = round(Int, div(fv-binedges[1], db))+1
                    if bnum < Nbins && bnum > 0
                        combined[bnum] += n
                    end
                end
            end

            nh = histograms[i,2]
            if nh > 0
                Nval = div(nh, 8) + zeros(Int, 8)
                for j=1:nh-sum(Nval)
                    Nval[rand(1:8)] += 1
                end
                values = linspace(binedges[i]+db/16., binedges[i+1]-db/16, 8)
                fvalues = f2(values)
                for (n, fv) in zip(Nval, fvalues)
                    bnum = round(Int, div(fv-binedges[1], db))+1
                    if bnum < Nbins && bnum > 0
                        combined[bnum] += n
                    end
                end
            end
        end
        return [f1,f2], combined
    end

    if N == 3
        curvesA, combinedA = matchspectra(histograms[:,1:2]; histrange=histrange)
        input = hcat(combinedA, histograms[:,3])
        curvesB, combined = matchspectra(input; histrange=histrange)
        f1 = compose_splinelog(curvesA[1], curvesB[1])
        f2 = compose_splinelog(curvesA[2], curvesB[1])
        f3 = curvesB[2]
        return [f1,f2,f3], combined
    end

    # Break the values vector up into equal-size groups (extra in N2 if N is odd)
    N1 = div(N,2)
    N2 = N - N1

    @assert N > 3
    histograms1 = histograms[:,1:N1]
    histograms2 = histograms[:,N1+1:end]

    curves1, combined1 = matchspectra(histograms1; histrange=histrange)
    curves2, combined2 = matchspectra(histograms2; histrange=histrange)

    @show N, typeof(combined1), eltype(combined1)
    @show N2, typeof(combined2), eltype(combined2)
    inputs = hcat(combined1, combined2)
    curves3, combined3 = matchspectra(inputs; histrange=histrange)
    finalcurves = []
    for c in curves1
        push!(finalcurves, compose_splinelog(c, curves3[1]))
    end
    for c in curves2
        push!(finalcurves, compose_splinelog(c, curves3[2]))
    end

    finalcurves, combined3
end

end # module


text:
using Test
using WordTokenizers

@testset "default behaviour" begin
	str = "The quick brown fox jumped over the lazy dog"
	tokenized = ["The", "quick", "brown", "fox", "jumped", "over", "the", "lazy", "dog"]

	@test tokenized == String.(rev_tokenize(str))
	@test str == rev_detokenize(tokenized)
end

@testset "multi-space" begin
    	str = "Some of   100,000 households (usually, a minority) ate breakfast.  "
    	tokenized = ["Some", "of", "100", "\ue302,\ue302", "000", "households", "(\ue302", "usually", "\ue302,", "a", "minority", "\ue302)", "ate", "breakfast", "\ue302."]

	@test tokenized == String.(rev_tokenize(str))
end

@testset "multi-byte character output" begin
	str = "Some of 100,000 households (usually, a minority) ate breakfast."
	tokenized = ["Some", "of", "100", "\ue302,\ue302", "000", "households", "(\ue302", "usually", "\ue302,", "a", "minority", "\ue302)", "ate", "breakfast", "\ue302."]

	@test tokenized == String.(rev_tokenize(str))
	@test str == rev_detokenize(tokenized)
end

@testset "multi-byte character input" begin
	str = "The quick brown fox jumped over the lazy dog ⌣⌣"
	tokenized = ["The", "quick", "brown", "fox", "jumped", "over", "the", "lazy", "dog", "⌣", "\ue302⌣"]

	@test tokenized == String.(rev_tokenize(str))
	@test str == rev_detokenize(tokenized)
end


text:
"""
    μ, σ = standardize!(X[, μ, σ; obsdim, operate_on])

or

    μ, σ = standardize!(D[, μ, σ; operate_on])

Standardize `X` along `obsdim` according to X = (X - μ) / σ.
If μ and σ are omitted they are computed such that variables have a mean of zero

`μ`         :  Vector or value describing the translation.
               Defaults to mean(X, 2)

`σ`         :  Vector or value describing the scale.
               Defaults to std(X, 2)

`obsdim`    :  Specify which axis corresponds to observations.
               Defaults to obsdim=2 (observations are columns of matrix)
               For DataFrames `obsdim` is obsolete and centering occurs
               column wise.

`operate_on`:  Specify the indices of columns or rows to be centered.
               Defaults to all columns/rows.
               For DataFrames this must be a vector of symbols, not indices
               E.g. `operate_on`=[1,3] will perform centering on columns
               with index 1 and 3 only (if obsdim=1, else rows 1 and 3)

Note on DataFrames:
Columns containing `NA` values are skipped.
Columns containing non numeric elements are skipped.

Examples:

    X = rand(4, 100)
    x = rand(10)
    D = DataFrame(A=rand(10), B=collect(1:10), C=[string(x) for x in 1:10])

    μ, σ = standardize!(X, obsdim=2)
    μ, σ = standardize!(X, ObsDim.First())
    μ, σ = standardize!(X, obsdim=1, operate_on=[1,3]
    μ, σ = standardize!(X, [7.0,8.0], [1,1], obsdim=1, operate_on=[1,3]
    μ, σ = standardize!(D)
    μ, σ = standardize!(D, operate_on=[:A,:B])
    μ, σ = standardize!(D, [-1,-1], [2,2], operate_on=[:A,:B])
"""
function standardize!(X, μ, σ; obsdim=LearnBase.default_obsdim(X), operate_on=default_scaleselection(X, convert(ObsDimension, obsdim)))
    standardize!(X, μ, σ, convert(ObsDimension, obsdim), operate_on)
end

function standardize!(X::AbstractArray{T,N}, μ, σ, ::ObsDim.Last, operate_on) where {T,N}
    standardize!(X, μ, σ, ObsDim.Constant{N}(), operate_on)
end

function standardize!(X; obsdim=LearnBase.default_obsdim(X), operate_on=default_scaleselection(X, convert(ObsDimension, obsdim)))
    standardize!(X, convert(ObsDimension, obsdim), operate_on)
end

function standardize!(X::AbstractArray{T,N}, ::ObsDim.Last, operate_on) where {T,N}
    standardize!(X, ObsDim.Constant{N}(), operate_on)
end

function standardize!(X::AbstractArray{T,N}, obsdim::ObsDim.Constant{M}, operate_on) where {T,N,M}
    μ = vec(mean(X, M))[operate_on]
    σ = vec(std(X, M))[operate_on]
    standardize!(X, μ, σ, obsdim, operate_on)
end

function standardize!(X::AbstractVector, ::ObsDim.Constant{M}, operate_on) where {M}
    μ = mean(X)
    σ = std(X)
    for i in operate_on
        X[i] = (X[i] - μ) / σ
    end
    μ, σ
end

function standardize!(X::AbstractMatrix, μ::AbstractVector, σ::AbstractVector, ::ObsDim.Constant{2}, operate_on)
    σ[σ .== 0] = 1
    nVars, nObs = size(X)
    for iObs in 1:nObs
        @inbounds for (i, iVar) in enumerate(operate_on)
            X[iVar, iObs] = (X[iVar, iObs] - μ[i]) / σ[i]
        end
    end
    μ, σ
end

function standardize!(X::AbstractMatrix, μ::AbstractVector, σ::AbstractVector, ::ObsDim.Constant{1}, operate_on)
    σ[σ .== 0] = 1
    nObs, nVars = size(X)
    for (i, iVar) in enumerate(operate_on)
        @inbounds for iObs in 1:nObs
            X[iObs, iVar] = (X[iObs, iVar] - μ[i]) / σ[i]
        end
    end
    μ, σ
end

function standardize!(X::AbstractVector, μ::AbstractVector, σ::AbstractVector, ::ObsDim.Constant{M}, operate_on) where {M}
    @inbounds for (i, iVar) in enumerate(operate_on)
        X[iVar] = (X[iVar] - μ[i]) / σ[i]
    end
    μ, σ
end

function standardize!(X::AbstractVector, μ::AbstractFloat, σ::AbstractFloat, ::ObsDim.Constant{M}, operate_on) where {M}
    @inbounds for i in 1:length(X)
        X[i] = (X[i] - μ) / σ
    end
    μ, σ
end

# --------------------------------------------------------------------
function standardize!(D::AbstractDataFrame; operate_on=default_scaleselection(D))
    standardize!(D, operate_on)
end

function standardize!(D::AbstractDataFrame, colnames::AbstractVector{Symbol})
    μ_vec = Float64[]
    σ_vec = Float64[]

    for colname in colnames
        if eltype(D[colname]) <: Real
            μ = mean(D[colname])
            σ = std(D[colname])
            if isna(μ)
                warn("Skipping \"$colname\" because it contains NA values")
                continue
            end
            standardize!(D, μ, σ, colname)
            push!(μ_vec, μ)
            push!(σ_vec, σ)
        else
            warn("Skipping \"$colname\" because data is not of type T <: Real.")
        end
    end
    μ_vec, σ_vec
end

function standardize!(D::AbstractDataFrame, μ::AbstractVector, σ::AbstractVector; operate_on=default_scaleselection(D))
    standardize!(D, μ, σ, operate_on)
end

function standardize!(D::AbstractDataFrame, μ::AbstractVector, σ::AbstractVector, colnames::AbstractVector{Symbol})
    for (icol, colname) in enumerate(colnames)
        standardize!(D, μ[icol], σ[icol], colname)
    end
    μ, σ
end

function standardize!(D::AbstractDataFrame, μ::Real, σ::Real, colname::Symbol)
    if any(isna, D[colname]) | !(eltype(D[colname]) <: Real)
        warn("Skipping \"$colname\" because it contains NA values or is not of type <: Real")
    else
        newcol::Vector{Float64} = convert(Vector{Float64}, D[colname])
        nobs = length(newcol)
        @inbounds for i in eachindex(newcol)
            newcol[i] = (newcol[i] - μ) / σ
       end
        D[colname] = newcol
    end
    μ, σ
end

"""
`StandardScaler` is used with the functions `fit()`, `transform()` and `fit_transform()`
to standardize data in a Matrix `X` or DataFrame according to Xnew = (X - μ) / σ.
After fitting a `StandardScaler` to one data set, it can be used to perform the same
transformation to a new set of data. E.g. fit the `StandardScaler` to your training
data and then apply the scaling to the test data at a later stage. (See examples below).

    fit(StandardScaler, X[, μ, σ; obsdim, operate_on])

    fit_transform(StandardScaler, X[, μ, σ; obsdim, operate_on])

`X`         :  Data of type Matrix or `DataFrame`.

`μ`         :  Vector or scalar describing the translation.
               Defaults to mean(X, obsdim)

`σ`         :  Vector or scalar describing the scale.
               Defaults to std(X, obsdim)

`obsdim`    :  Specify which axis corresponds to observations.
               Defaults to obsdim=2 (observations are columns of matrix)
               For DataFrames `obsdim` is obsolete and rescaling occurs
               column wise.

`operate_on`:  Specify the indices of columns or rows to be centered.
               Defaults to all columns/rows.
               For DataFrames this must be a vector of symbols, not indices.
               E.g. `operate_on`=[1,3] will perform centering on columns
               with index 1 and 3 only (if obsdim=1, else rows 1 and 3)

Note on DataFrames:
Columns containing `NA` values are skipped.
Columns containing non numeric elements are skipped.

Examples:

    Xtrain = rand(100, 4)
    Xtest  = rand(10, 4)
    x = rand(4)
    Dtrain = DataFrame(A=rand(10), B=collect(1:10), C=[string(x) for x in 1:10])
    Dtest = DataFrame(A=rand(10), B=collect(1:10), C=[string(x) for x in 1:10])

    scaler = fit(StandardScaler, Xtrain)
    scaler = fit(StandardScaler, Xtrain, obsdim=1)
    scaler = fit(StandardScaler, Xtrain, obsdim=1, operate_on=[1,3])
    transform(Xtest, scaler)
    transform!(Xtest, scaler)
    transform(x, scaler)
    transform!(x, scaler)

    scaler = fit(StandardScaler, Dtrain)
    scaler = fit(StandardScaler, Dtrain, operate_on=[:A,:B])
    transform(Dtest, scaler)
    transform!(Dtest, scaler)

    Xscaled, scaler = fit_transform(StandardScaler, X, obsdim=1, operate_on=[1,2,4])
    scaler = fit_transform!(StandardScaler, X, obsdim=1, operate_on=[1,2,4])

Note that for `transform!` the data matrix `X` has to be of type <: AbstractFloat
as the scaling occurs inplace. (E.g. cannot be of type Matrix{Int64}). This is not
the case for `transform` however.
For `DataFrames` `transform!` can be used on columns of type <: Integer.
"""
struct StandardScaler{T<:Real,U<:Real,I,M}
    offset::Vector{T}
    scale::Vector{U}
    obsdim::ObsDim.Constant{M}
    operate_on::Vector{I}
end

function StandardScaler(X::AbstractArray{T,M}; obsdim=LearnBase.default_obsdim(X), operate_on=default_scaleselection(X, convert(ObsDimension, obsdim))) where {T<:Real,M}
    StandardScaler(X, convert(ObsDimension, obsdim), operate_on)
end

function StandardScaler(X::AbstractArray{T,M}, ::ObsDim.Last, operate_on) where {T<:Real,M}
    StandardScaler(X, ObsDim.Constant{M}(), operate_on)
end

function StandardScaler(X::AbstractArray{T,N}, obsdim::ObsDim.Constant{M}, operate_on::AbstractVector) where {T<:Real,N,M}
    offset = vec(mean(X,M))[operate_on]
    scale = vec(std(X, M))[operate_on]
    StandardScaler(offset, scale, obsdim, operate_on)
end

function StandardScaler(D::AbstractDataFrame; operate_on=default_scaleselection(D))
    StandardScaler(D, operate_on)
end

function StandardScaler(D::AbstractDataFrame, operate_on::Vector{Symbol})
    colnames = valid_columns(D, operate_on)
    offset = Float64[mean(D[colname]) for colname in colnames]
    scale = Float64[std(D[colname]) for colname in colnames]
    StandardScaler(offset, scale, ObsDim.Constant{1}(), colnames)
end

function StandardScaler(D::AbstractDataFrame, offset, scale; operate_on=default_scaleselection(D))
    colnames = valid_columns(D)
    StandardScaler(offset, scale, ObsDim.Constant{1}(), colnames)
end

function StatsBase.fit(::Type{StandardScaler}, X::AbstractMatrix{T}; obsdim=LearnBase.default_obsdim(X), operate_on=default_scaleselection(X, convert(ObsDimension, obsdim))) where {T<:Real}
    StandardScaler(X, convert(ObsDimension, obsdim), operate_on)
end

function fit_transform(::Type{StandardScaler}, X::AbstractMatrix{T}; obsdim=LearnBase.default_obsdim(X), operate_on=default_scaleselection(X, convert(ObsDimension, obsdim))) where {T<:Real}
    scaler = StandardScaler(X, convert(ObsDimension, obsdim), operate_on)
    Xnew = transform(X, scaler)
    return Xnew, scaler
end

function fit_transform!(::Type{StandardScaler}, X::AbstractMatrix{T}; obsdim=LearnBase.default_obsdim(X), operate_on=default_scaleselection(X, convert(ObsDimension, obsdim))) where {T<:Real}
    scaler = StandardScaler(X, convert(ObsDimension, obsdim), operate_on)
    transform!(X, scaler)
    return scaler
end

function StatsBase.fit(::Type{StandardScaler}, D::AbstractDataFrame; operate_on=default_scaleselection(D))
    StandardScaler(D, operate_on)
end

function fit_transform(::Type{StandardScaler}, D::AbstractDataFrame; operate_on=default_scaleselection(D))
    scaler = StandardScaler(D, operate_on)
    Dnew = transform(D, scaler)
    return Dnew, scaler
end

function fit_transform!(::Type{StandardScaler}, D::AbstractDataFrame; operate_on=default_scaleselection(D))
    scaler = StandardScaler(D, operate_on)
    transform!(D, scaler)
    return scaler
end

function transform!(X::AbstractArray{T,N}, cs::StandardScaler) where {T<:AbstractFloat,N}
    standardize!(X, cs.offset, cs.scale, cs.obsdim, cs.operate_on)
    X
end

function transform!(D::AbstractDataFrame, cs::StandardScaler)
    standardize!(D, cs.offset, cs.scale, cs.operate_on)
    D
end

function transform(X::AbstractArray{T,N}, cs::StandardScaler) where {T<:AbstractFloat,N}
    Xnew = deepcopy(X)
    transform!(Xnew, cs)
end

function transform(X::AbstractArray{T,N}, cs::StandardScaler) where {T<:Real,N}
    Xnew = convert(AbstractArray{Float64, N}, X)
    transform!(Xnew, cs)
    Xnew
end

function transform(D::AbstractDataFrame, cs::StandardScaler)
    Dnew = deepcopy(D)
    transform!(Dnew, cs)
    Dnew
end


text:
"""
$(TYPEDEF)

When defining problem-specific cache, make it a subtype of this.
"""
abstract type AbstractExtraILMCache end

"""
$(TYPEDEF)

A cache of operators and storage data for use in surface operations. Constructed
with [`SurfaceScalarCache`](@ref) or [`SurfaceVectorCache`](@ref).
"""
struct BasicILMCache{N,SCA<:AbstractScalingType,GCT,ND,BLT<:BodyList,NT<:VectorData,
                      DST<:ScalarData,REGT<:Regularize,
                      RSNT<:RegularizationMatrix,ESNT<:InterpolationMatrix,
                      RT<:RegularizationMatrix,ET<:InterpolationMatrix,
                      RCT, ECT, RDT, EDT,
                      LT<:CartesianGrids.Laplacian,GVT,GNT,GDT,SVT,SDT,SST} <: AbstractBasicCache{N,GCT}

    # Grid
    g :: PhysicalGrid{ND}

    # Bodies
    bl :: BLT

    # Normals
    nrm :: NT

    # Areas
    ds :: DST

    # Regularization operator
    regop :: REGT

    # Regularization and interpolation of tensor product data (gsnorm/snorm)
    Rsn :: RSNT
    Esn :: ESNT

    # Regularization and interpolation of basic data (gdata/sdata)
    R :: RT
    E :: ET

    # Regularization and interpolation of grid curl data (gcurl/sdata)
    Rcurl :: RCT
    Ecurl :: ECT

    # Regularization and interpolation of grid div data (gdiv/sdata)
    Rdiv :: RDT
    Ediv :: EDT

    # Laplacian (with no coefficient)
    L :: LT

    # For holding the grid data comprised of tensor product of basic grid data
    # and regularized normals (e.g., scalar -> Edges{Primal}, vector -> EdgeGradient)
    gsnorm_cache :: GVT
    gsnorm2_cache :: GVT

    # For holding the curl of the basic grid data (always Nodes{Dual})
    gcurl_cache :: GNT

    # For holding the div of the basic grid data (always Nodes{Primal})
    gdiv_cache :: GDT

    # For holding the basic data type (e.g., scalar -> Nodes{Primal}, vector -> Edges{Primal})
    gdata_cache :: GCT

    # For holding the surface data comprised of tensor product of basic data
    # and normals (e.g, scalar -> VectorData, vector -> TensorData)
    snorm_cache :: SVT
    snorm2_cache :: SVT

    # For holding the basic surface data (e.g., scalar -> ScalarData, vector -> VectorData)
    sdata_cache :: SDT

    # For holding surface scalar data (always ScalarData)
    sscalar_cache :: SST
end

RigidBodyTools.numpts(::BasicILMCache{N}) where {N} = N
RigidBodyTools.numpts(::AbstractBasicCache{N}) where {N} = N
scalingtype(::BasicILMCache{N,SCA}) where {N,SCA} = SCA
cache_datatype(::BasicILMCache{N,SCA,GCT}) where {N,SCA,GCT} = GCT <: Nodes ? :scalar : :vector
cache_datatype(::AbstractBasicCache{N,GCT}) where {N,GCT} = GCT <: Nodes ? :scalar : :vector
type_curlreg(::BasicILMCache{N,SCA,GCT,ND,BLT,NT,DST,REGT,RSNT,ESNT,RT,ET,RCT,ECT,RDT,EDT,LT,GVT,GNT,GDT,SVT,SDT,SST}) where {N,SCA,GCT,ND,BLT,NT,DST,REGT,RSNT,ESNT,RT,ET,RCT,ECT,RDT,EDT,LT,GVT,GNT,GDT,SVT,SDT,SST} = RCT
type_divreg(::BasicILMCache{N,SCA,GCT,ND,BLT,NT,DST,REGT,RSNT,ESNT,RT,ET,RCT,ECT,RDT,EDT,LT,GVT,GNT,GDT,SVT,SDT,SST}) where {N,SCA,GCT,ND,BLT,NT,DST,REGT,RSNT,ESNT,RT,ET,RCT,ECT,RDT,EDT,LT,GVT,GNT,GDT,SVT,SDT,SST} = RDT
type_sdata(::BasicILMCache{N,SCA,GCT,ND,BLT,NT,DST,REGT,RSNT,ESNT,RT,ET,RCT,ECT,RDT,EDT,LT,GVT,GNT,GDT,SVT,SDT,SST}) where {N,SCA,GCT,ND,BLT,NT,DST,REGT,RSNT,ESNT,RT,ET,RCT,ECT,RDT,EDT,LT,GVT,GNT,GDT,SVT,SDT,SST} = SDT
type_sscalar(::BasicILMCache{N,SCA,GCT,ND,BLT,NT,DST,REGT,RSNT,ESNT,RT,ET,RCT,ECT,RDT,EDT,LT,GVT,GNT,GDT,SVT,SDT,SST}) where {N,SCA,GCT,ND,BLT,NT,DST,REGT,RSNT,ESNT,RT,ET,RCT,ECT,RDT,EDT,LT,GVT,GNT,GDT,SVT,SDT,SST} = SST

"""
    SurfaceScalarCache(g::PhysicalGrid[,scaling=IndexScaling])

Create a cache of type `BasicILMCache` with scalar grid data, using the grid specified
in `g`, with no immersed points. The keyword `scaling`
can be used to set the scaling in the operations.
By default, it is set to `IndexScaling` which sets the differential operators
to be only differencing operators. By using `scaling = GridScaling`, then the grid and
 spacings are accounted for and differential operators are scaled by this spacing.
 The keyword `phys_params` can be used to supply physical parameters.
""" SurfaceScalarCache(::PhysicalGrid)

"""
    SurfaceVectorCache(g::PhysicalGrid[,scaling=IndexScaling])

Create a cache of type `BasicILMCache` with vector grid data, with no immersed points.
See [`SurfaceScalarCache`](@ref) for details.
""" SurfaceVectorCache(::PhysicalGrid)

"""
    SurfaceVectorCache(body::Body/BodyList,g::PhysicalGrid[,ddftype=CartesianGrids.Yang3][,scaling=Gridcaling])

Create a cache of type `BasicILMCache`, holding operators and storage data
for use in immersed layer operations on vector data. See [`SurfaceScalarCache`](@ref)
for details.
""" SurfaceVectorCache(::Union{Body,BodyList},::PhysicalGrid)

"""
    SurfaceScalarCache(body::Body/BodyList,g::PhysicalGrid[,ddftype=CartesianGrids.Yang3][,scaling=GridScaling])

Create a cache of type `BasicILMCache`, holding operators and storage data
for use in immersed layer operations on scalar data. This is sometimes called from within`ILMSystem` rather than directly.

The `body` can be of type `Body` or `BodyList`. The keyword `scaling` can be used to set the scaling in the operations.
By default, it is set to `IndexScaling` which sets the regularization and interpolation to
be symmetric matrices (i.e., interpolation is the adjoint of regularization with
  respect to a vector dot product), and the vector calculus operations on the grid
  are simple differences. By using `scaling = GridScaling`, then the grid and
  point spacings are accounted for. Interpolation and regularization are adjoints
  with respect to inner products based on discretized surface and volume integrals,
  and vector calculus operations are scaled by the grid spacing.
""" SurfaceScalarCache(::Union{Body,BodyList},::PhysicalGrid)

"""
    SurfaceScalarCache(X::VectorData,g::PhysicalGrid[,ddftype=CartesianGrids.Yang3][,scaling=GridScaling])

Create a cache of type `BasicILMCache`, holding operators and storage data
for use in immersed layer operations on scalar data. The `X` specifies the
is assumed to hold the endpoints of the immersed surface segments, and `g` the physical grid.
""" SurfaceScalarCache(::VectorData,::ScalarData,::VectorData,::PhysicalGrid)

"""
    SurfaceVectorCache(X::VectorData,g::PhysicalGrid[,ddftype=CartesianGrids.Yang3][,scaling=GridScaling])

Create a cache of type `BasicILMCache`, holding operators and storage data
for use in immersed layer operations on vector data. See [`SurfaceScalarCache`](@ref)
for details.
""" SurfaceVectorCache(::VectorData,::ScalarData,::VectorData,::PhysicalGrid)

function SurfaceScalarCache(bl::BodyList,a::ScalarData{N},nrm::VectorData{N},g::PhysicalGrid;
                              ddftype = DEFAULT_DDF,
                              scaling = DEFAULT_SCALING,
                              dtype = DEFAULT_DATA_TYPE) where {N}

	X = points(bl)
  sscalar_cache = nothing
	sdata_cache = ScalarData(X, dtype = dtype)
	snorm_cache = VectorData(X, dtype = dtype)

	gsnorm_cache = Edges(Primal,size(g), dtype = dtype)
	gdata_cache = Nodes(Primal,size(g), dtype = dtype)
  gcurl_cache = Nodes(Dual,size(g), dtype = dtype)
  gdiv_cache = nothing #Nodes(Primal,size(g), dtype = dtype)

	_surfacecache(bl,X,a,nrm,g,ddftype,scaling,sdata_cache,snorm_cache,sscalar_cache,gsnorm_cache,gcurl_cache,gdiv_cache,gdata_cache;dtype=dtype)

end

function SurfaceVectorCache(bl::BodyList,a::ScalarData{N},nrm::VectorData{N},g::PhysicalGrid;
                              ddftype = DEFAULT_DDF,
                              scaling = DEFAULT_SCALING,
                              dtype = DEFAULT_DATA_TYPE) where {N}

	X = points(bl)
  sscalar_cache = ScalarData(X, dtype = dtype)
	sdata_cache = VectorData(X, dtype = dtype)
	snorm_cache = TensorData(X, dtype = dtype)

	gsnorm_cache = EdgeGradient(Primal,size(g), dtype = dtype)
  gdata_cache = Edges(Primal,size(g), dtype = dtype)
	gcurl_cache = Nodes(Dual,size(g), dtype = dtype)
  gdiv_cache = Nodes(Primal,size(g), dtype = dtype)

	_surfacecache(bl,X,a,nrm,g,ddftype,scaling,sdata_cache,snorm_cache,sscalar_cache,gsnorm_cache,gcurl_cache,gdiv_cache,gdata_cache;dtype=dtype)

end

for f in [:SurfaceScalarCache, :SurfaceVectorCache]
  @eval $f(body::Body,g::PhysicalGrid; kwargs...) =
        $f(BodyList([body]),areas(body),normals(body),g; kwargs...)

  @eval $f(bl::BodyList,g::PhysicalGrid; kwargs...) =
        $f(bl,areas(bl),normals(bl),g; kwargs...)

  @eval $f(g::PhysicalGrid;kwargs...) =
        $f(BodyList(),ScalarData(0),VectorData(0),g; kwargs...)

  @eval function $f(X::VectorData,g::PhysicalGrid; kwargs...)
          x = Vector{Float64}(undef,length(X.u))
          y = Vector{Float64}(undef,length(X.v))
          x .= X.u
          y .= X.v
          $f(BasicBody(x,y),g; kwargs...)
  end

end

function Base.show(io::IO, H::BasicILMCache{N,SCA}) where {N,SCA}
    println(io, "Surface cache with scaling of type $SCA")
    println(io, "  $N point data of type $(typeof(H.sdata_cache))")
    println(io, "  Grid data of type $(typeof(H.gdata_cache))")
end

function _surfacecache(bl::BodyList,X::VectorData{N},a,nrm,g::PhysicalGrid{ND},ddftype,scaling,
                      sdata_cache,snorm_cache,sscalar_cache,gsnorm_cache,gcurl_cache,gdiv_cache,gdata_cache;dtype=Float64) where {N,ND}

  regop = _get_regularization(X,a,g,ddftype,scaling)
  Rsn = _regularization_matrix(regop,snorm_cache,gsnorm_cache)
  Esn = _interpolation_matrix(regop, gsnorm_cache, snorm_cache)

  R = _regularization_matrix(regop,sdata_cache,gdata_cache )
  E = _interpolation_matrix(regop, gdata_cache,sdata_cache)

  Rcurl = _regularization_matrix(regop,sscalar_cache,gcurl_cache )
  Ecurl = _interpolation_matrix(regop, gcurl_cache,sscalar_cache)

  Rdiv = _regularization_matrix(regop,sscalar_cache,gdiv_cache )
  Ediv = _interpolation_matrix(regop, gdiv_cache,sscalar_cache)

  coeff_factor = 1.0
  with_inverse = true
  L = _get_laplacian(gdata_cache,coeff_factor,g,with_inverse,scaling;dtype=dtype)

  return BasicILMCache{N,scaling,typeof(gdata_cache),ND,typeof(bl),typeof(nrm),typeof(a),typeof(regop),
                       typeof(Rsn),typeof(Esn),typeof(R),typeof(E),typeof(Rcurl),typeof(Ecurl),typeof(Rdiv),typeof(Ediv),typeof(L),
                       typeof(gsnorm_cache),typeof(gcurl_cache),typeof(gdiv_cache),typeof(snorm_cache),typeof(sdata_cache),typeof(sscalar_cache)}(
                       g,bl,nrm,a,regop,Rsn,Esn,R,E,Rcurl,Ecurl,Rdiv,Ediv,L,
                       _similar(gsnorm_cache),_similar(gsnorm_cache),
                       _similar(gcurl_cache),_similar(gdiv_cache),_similar(gdata_cache),
                       _similar(snorm_cache),_similar(snorm_cache),_similar(sdata_cache),_similar(sscalar_cache))

end

#=
Point collection cache
=#

struct PointCollectionCache{N,GCT,ND,PT,REGT<:Regularize,SST} <: AbstractBasicCache{N,GCT}
    pts :: PT
    g :: PhysicalGrid{ND}
    regop :: REGT
    gdata_cache :: GCT
    sdata_cache :: SST
end

function ScalarPointCollectionCache(X::VectorData{N},g::PhysicalGrid{ND};ddftype = DEFAULT_DDF,dtype = DEFAULT_DATA_TYPE) where {N,ND}
    gdata_cache = Nodes(Primal,size(g),dtype = dtype)
    sdata_cache = ScalarData(X,dtype = dtype)
    return _pointcollectioncache(X,g,gdata_cache,sdata_cache,ddftype)
end

function VectorPointCollectionCache(X::VectorData{N},g::PhysicalGrid{ND};ddftype = DEFAULT_DDF,dtype = DEFAULT_DATA_TYPE) where {N,ND}
    gdata_cache = Edges(Primal,size(g),dtype = dtype)
    sdata_cache = VectorData(X,dtype = dtype)
    return _pointcollectioncache(X,g,gdata_cache,sdata_cache,ddftype)
end

function _pointcollectioncache(X::VectorData{N},g::PhysicalGrid{ND},gdata_cache::GCT,sdata_cache::SDT,ddftype) where {N,ND,GCT,SDT}
  regop = _get_regularization(X,g,ddftype)
  return PointCollectionCache{N,GCT,ND,typeof(X),typeof(regop),SDT}(X,g,regop,gdata_cache,sdata_cache)
end

#=
Convenience functions
=#

@inline CartesianGrids.cellsize(s::AbstractBasicCache) = cellsize(s.g)
@inline Base.length(s::AbstractBasicCache{N}) where {N} = N
Base.eltype(s::AbstractBasicCache) = eltype(s.gdata_cache)

# Standardize the regularization
for ddf in [:Yang3,:Roma,:Goza,:Witchhat,:M3,:M4prime]
  @eval _get_regularization(X::VectorData{N},a::ScalarData{N},g::PhysicalGrid,::Type{CartesianGrids.$ddf},::Type{GridScaling};filter=false) where {N} =
     Regularize(X,cellsize(g),I0=origin(g),weights=a.data,ddftype=CartesianGrids.$ddf,filter=filter)

  @eval _get_regularization(X::VectorData{N},a::ScalarData{N},g::PhysicalGrid,::Type{CartesianGrids.$ddf},::Type{IndexScaling};filter=false) where {N} =
     Regularize(X,cellsize(g),I0=origin(g),issymmetric=true,ddftype=CartesianGrids.$ddf,filter=filter)

  @eval _get_regularization(X::VectorData{N},g::PhysicalGrid,::Type{CartesianGrids.$ddf}) where {N} =
     Regularize(X,cellsize(g),I0=origin(g),ddftype=CartesianGrids.$ddf,filter=false)
end

_get_regularization(body::Union{Body,BodyList},args...;kwargs...) = _get_regularization(VectorData(collect(body)),areas(body),args...;kwargs...)

# Standardize the Laplacian
_get_laplacian(a,coeff_factor::Real,g::PhysicalGrid,with_inverse,::Type{IndexScaling};dtype=Float64) =
               CartesianGrids.plan_laplacian(a,with_inverse=with_inverse,factor=coeff_factor,dtype=dtype)
_get_laplacian(a,coeff_factor::Real,g::PhysicalGrid,with_inverse,::Type{GridScaling};dtype=Float64) =
               CartesianGrids.plan_laplacian(a,with_inverse=with_inverse,factor=coeff_factor/cellsize(g)^2,dtype=dtype)

# This is needed to stabilize the type-unstable `RegularizationMatrix` function in
# CartesianGrids
function _regularization_matrix(regop::Regularize,src::PointData,trg::GridData)
    if regop._issymmetric
      R, _ = RegularizationMatrix(regop, src, trg)
    else
      R = RegularizationMatrix(regop, src, trg)
    end
    return R
end

@inline _regularization_matrix(regop::Regularize,src::PointData,::Nothing) = nothing

@inline _regularization_matrix(regop::Regularize,::Nothing,trg) = nothing

@inline _interpolation_matrix(regop::Regularize,src::GridData,trg::PointData) =
        InterpolationMatrix(regop,src,trg)

@inline _interpolation_matrix(regop::Regularize,::Nothing,trg::PointData) = nothing

@inline _interpolation_matrix(regop::Regularize,src,::Nothing) = nothing

# APIs to generate regularization, interpolation matrices and Laplacians not generated
# in the basic cache
"""
    RegularizationMatrix(cache::BasicILMCache,src::PointData,trg::GridData)

Create a regularization matrix for regularizing point data of type `src` to
grid data of type `trg`. (Both `src` and `trg` must be appropriately sized
for the grid and points in `cache`.)
"""
CartesianGrids.RegularizationMatrix(cache::BasicILMCache,src::PointData,trg::GridData) =
    _regularization_matrix(cache.regop,src,trg)

"""
    InterpolationMatrix(cache::BasicILMCache,src::GridData,trg::PointData)

Create a interpolation matrix for regularizing grid data of type `src` to
point data of type `trg`. (Both `src` and `trg` must be appropriately sized
for the grid and points in `cache`.)
"""
CartesianGrids.InterpolationMatrix(cache::BasicILMCache,src::GridData,trg::PointData) =
    _interpolation_matrix(cache.regop,src,trg)

"""
    Laplacian(cache::BasicILMCache,src::GridData,coeff_factor::Real[,with_inverse=true])

Create an invertible Laplacian operator for data `src` (of the same size as the grid in `cache`),
using the index or grid scaling associated with `cache`. The operator is pre-multiplied
by the factor `coeff_factor`.
"""
Laplacian(cache::BasicILMCache{N,SCA},src::GridData,coeff_factor; with_inverse=true, dtype=Float64) where {N,SCA} =
    _get_laplacian(src,coeff_factor,cache.g,with_inverse,SCA;dtype=dtype)

# Some utilities to get the DDF type of the cache
_ddf_type(::DDF{DT}) where {DT} = DT
_ddf_type(R::Regularize) = _ddf_type(R.ddf)
_ddf_type(cache::BasicILMCache) = _ddf_type(cache.regop)

# Getting the list of first indices for each body, for partitioning of surface vectors
_firstindices(b::Body) = [1,length(body)+1]
_firstindices(bl::BodyList) = [map(i -> first(getrange(bl,i)),1:length(bl)); numpts(bl)+1]

## Obtaining copies of the grid and surface data
"""
    similar_grid(::BasicILMCache)

Get a `similar` copy of the basic grid data in the cache.
"""
@inline similar_grid(cache::AbstractBasicCache,kwargs...) = _similar(cache.gdata_cache,kwargs...)

"""
    similar_gridgrad(::BasicILMCache)

Get a `similar` copy of the gradient of the grid data in the cache.
"""
@inline similar_gridgrad(cache::BasicILMCache,kwargs...) = _similar(cache.gsnorm_cache,kwargs...)

"""
    similar_gridcurl(::BasicILMCache)

Get a `similar` copy of the grid curl field data in the cache.
"""
@inline similar_gridcurl(cache::BasicILMCache,kwargs...) = _similar(cache.gcurl_cache,kwargs...)

"""
    similar_griddiv(::BasicILMCache)

Get a `similar` copy of the grid div field data in the cache.
"""
@inline similar_griddiv(cache::BasicILMCache,kwargs...) = _similar(cache.gdiv_cache,kwargs...)

"""
    similar_gridgradcurl(::BasicILMCache)

Get a `similar` copy of the grid gradient-of-curl field data in the cache.
"""
@inline similar_gridgradcurl(cache::BasicILMCache;element_type=eltype(cache)) = Edges(Dual,cache.gdata_cache,dtype=element_type)

"""
    similar_surface(::BasicILMCache)

Get a `similar` copy of the basic surface point data in the cache.
"""
@inline similar_surface(cache::AbstractBasicCache,kwargs...) = _similar(cache.sdata_cache,kwargs...)

"""
    similar_surfacescalar(::BasicILMCache)

Get a `similar` copy of the surface scalar point data in the cache. This
is only used for vector-type caches. Otherwise, [`similar_surface`](@ref) should
be used.
"""
@inline similar_surfacescalar(cache::AbstractBasicCache,kwargs...) = _similar(cache.sscalar_cache,kwargs...)

_similar(::Nothing;kwargs...) = nothing
_similar(a;kwargs...) = similar(a;kwargs...)

"""
    zeros_grid(::BasicILMCache)

Get an instance of the basic grid data in the cache, with values set to zero.
"""
@inline zeros_grid(cache::AbstractBasicCache,kwargs...) = _zero(cache.gdata_cache,kwargs...)

"""
    zeros_gridgrad(::BasicILMCache,dim)

Get an instance of the gradient of the grid data in the cache, , in direction `dim`,
with values set to zero. If the data are of type `TensorGridData`, then
`dim` takes values from 1 to 2^2.
"""
@inline zeros_gridgrad(cache::BasicILMCache,kwargs...) = _zero(cache.gsnorm_cache,kwargs...)

"""
    zeros_gridcurl(::BasicILMCache)

Get an instance of the grid curl field data in the cache, with values set to zero.
"""
@inline zeros_gridcurl(cache::BasicILMCache,kwargs...) = _zero(cache.gcurl_cache,kwargs...)

"""
    zeros_griddiv(::BasicILMCache)

Get an instance of the grid div field data in the cache, with values set to zero.
"""
@inline zeros_griddiv(cache::BasicILMCache,kwargs...) = _zero(cache.gdiv_cache,kwargs...)

"""
    zeros_gridgradcurl(::BasicILMCache)

Get an instance of the grid gradient-of-curl field data in the cache, with values set to zero.
"""
@inline zeros_gridgradcurl(cache::BasicILMCache;element_type=eltype(cache)) = Edges(Dual,cache.gdata_cache,dtype=element_type)

"""
    zeros_surface(::BasicILMCache)

Get an instance of the basic surface point data in the cache, with values set to zero.
"""
@inline zeros_surface(cache::AbstractBasicCache,kwargs...) = _zero(cache.sdata_cache,kwargs...)

"""
    zeros_surfacescalar(::BasicILMCache)

Get an instance of the surface scalar point data in the cache, with values set to zero.
This is only used for vector-type caches. Otherwise, [`zeros_surface`](@ref) should
be used.
"""
@inline zeros_surfacescalar(cache::AbstractBasicCache,kwargs...) = _zero(cache.sscalar_cache,kwargs...)

_zero(a;kwargs...) = zero(a;kwargs...)
_zero(::Nothing;kwargs...) = nothing

"""
    ones_grid(::BasicILMCache)

Get an instance of the basic grid data in the cache, with values set to unity.
"""
@inline ones_grid(cache::AbstractBasicCache;kwargs...) = ones(cache.gdata_cache;kwargs...)

"""
    ones_grid(::BasicILMCache,dim)

For a vector-type cache, get an instance of the basic grid data in the cache, with values set to unity
in dimension `dim`.
"""
@inline ones_grid(cache::AbstractBasicCache{N,GCT},dim;kwargs...) where {N,GCT<:Edges} = ones(cache.gdata_cache,dim;kwargs...)

"""
    ones_gridgrad(::BasicILMCache,dim)

Get an instance of the gradient of the grid data in the cache, in direction `dim`,
with values set to unity. If the data are of type `TensorGridData`, then
`dim` takes values from 1 to 2^2.
"""
@inline ones_gridgrad(cache::BasicILMCache,dim;kwargs...) = ones(cache.gsnorm_cache,dim;kwargs...)

"""
    ones_gridcurl(::BasicILMCache)

Get an instance of the grid curl field data in the cache, with values set to unity.
"""
@inline ones_gridcurl(cache::BasicILMCache,kwargs...) = _ones(cache.gcurl_cache,kwargs...)

"""
    ones_griddiv(::BasicILMCache)

Get an instance of the grid div field data in the cache, with values set to unity.
"""
@inline ones_griddiv(cache::BasicILMCache,kwargs...) = _ones(cache.gdiv_cache,kwargs...)

"""
    ones_gridgradcurl(::BasicILMCache)

Get an instance of the grid gradient-of-curl field data in the cache, with values set to unity.
"""
@inline ones_gridgradcurl(cache::BasicILMCache;element_type=eltype(cache)) =
      (d = Edges(Dual,cache.gdata_cache,dtype=element_type); fill!(d,one(element_type)))

"""
    ones_surface(::BasicILMCache)

Get an instance of the basic surface point data in the cache, with values set to unity.
"""
@inline ones_surface(cache::AbstractBasicCache;kwargs...) = ones(cache.sdata_cache;kwargs...)

"""
    ones_surface(::BasicILMCache,dim)

Get an instance of the basic surface point data in the cache, with values set to unity
in dimension `dim`. This only works for a vector-type cache.
"""
@inline ones_surface(cache::AbstractBasicCache{N,GCT},dim;kwargs...) where {N,GCT<:Edges} = ones(cache.sdata_cache,dim;kwargs...)

"""
    ones_surfacescalar(::BasicILMCache)

Get an instance of the surface scalar point data in the cache, with values set to unity.
This is only used for vector-type caches. Otherwise, [`ones_surface`](@ref) should
be used.
"""
@inline ones_surfacescalar(cache::AbstractBasicCache;kwargs...) = _ones(cache.sscalar_cache;kwargs...)

_ones(a...;kwargs...) = ones(a...;kwargs...)
_ones(::Nothing;kwargs...) = nothing

"""
    x_grid(::BasicILMCache)

Return basic grid data filled with the grid `x` coordinate.
If the grid data is scalar, then the output of this function is
of the same type. If the grid
data is vector, then the output is of type `Edges{Primal}`,
and the coordinates of each component are in `u`, `v` fields.
"""
function x_grid(cache::BasicILMCache{N,SCA,GT}) where {N,SCA,GT<:Nodes{Primal}}
    xc, _ = coordinates(cache.gdata_cache,cache.g)
    p = zeros_grid(cache)
    p .= xc
    return p
end

function x_grid(cache::BasicILMCache{N,SCA,GT}) where {N,SCA,GT<:Edges{Primal}}
    p = zeros_grid(cache)
    xu, _ = coordinates(p.u,cache.g)
    xv, _ = coordinates(p.v,cache.g)
    p.u .= xu
    p.v .= xv
    return p
end

"""
    x_gridcurl(::BasicILMCache)

Return basic grid curl field data filled with the grid `x` coordinate
"""
function x_gridcurl(cache::BasicILMCache)
    xc, _ = coordinates(cache.gcurl_cache,cache.g)
    p = zeros_gridcurl(cache)
    p .= xc
end

"""
    x_griddiv(::BasicILMCache)

Return basic grid div field data filled with the grid `x` coordinate
"""
@inline x_griddiv(cache::BasicILMCache) = _x_griddiv(cache,Val(cache_datatype(cache)))

_x_griddiv(cache,::Val{:scalar}) = nothing

function _x_griddiv(cache,::Val{:vector})
    xc, _ = coordinates(cache.gdiv_cache,cache.g)
    p = zeros_griddiv(cache)
    p .= xc
end

"""
    x_gridgrad(::BasicILMCache)

Return basic grid gradient field data filled with the grid `x` coordinate.
If the grid data is scalar, then the output of this function is
of `Edges{Primal}` type and the coordinate field for each component
can be accessed with the `u` and `v` field, respectively. If the grid
data is vector, then the output is of type `EdgeGradient{Primal}`,
and the coordinates are in `dudx`, `dudy`, `dvdx`, and `dvdy` fields.
"""
function x_gridgrad(cache::BasicILMCache{N,SCA,GT}) where {N,SCA,GT<:Nodes{Primal}}
    p = zeros_gridgrad(cache)
    xu, _ = coordinates(p.u,cache.g)
    xv, _ = coordinates(p.v,cache.g)
    p.u .= xu
    p.v .= xv
    return p
end

function x_gridgrad(cache::BasicILMCache{N,SCA,GT}) where {N,SCA,GT<:Edges{Primal}}
    p = zeros_gridgrad(cache)
    xdudx, _ = coordinates(p.dudx,cache.g)
    xdudy, _ = coordinates(p.dudy,cache.g)
    p.dudx .= xdudx
    p.dudy .= xdudy
    p.dvdx .= xdudy
    p.dvdy .= xdudx
    return p
end

"""
    y_grid(::BasicILMCache)

Return basic grid data filled with the grid `y` coordinate.
If the grid data is scalar, then the output of this function is
of the same type. If the grid
data is vector, then the output is of type `Edges{Primal}`,
and the coordinates of each component are in `u`, `v` fields.
"""
function y_grid(cache::BasicILMCache{N,SCA,GT}) where {N,SCA,GT<:Nodes{Primal}}
    _, yc  = coordinates(cache.gdata_cache,cache.g)
    p = zeros_grid(cache)
    p .= yc
    return p
end

function y_grid(cache::BasicILMCache{N,SCA,GT}) where {N,SCA,GT<:Edges{Primal}}
    p = zeros_grid(cache)
    _, yu  = coordinates(p.u,cache.g)
    _, yv  = coordinates(p.v,cache.g)
    p.u .= yu'
    p.v .= yv'
    return p
end

"""
    y_gridcurl(::BasicILMCache)

Return basic grid curl field data filled with the grid `y` coordinate
"""
function y_gridcurl(cache::BasicILMCache)
    _,yc = coordinates(cache.gcurl_cache,cache.g)
    p = zeros_gridcurl(cache)
    p .= yc'
end

"""
    y_griddiv(::BasicILMCache)

Return basic grid div field data filled with the grid `y` coordinate
"""
@inline y_griddiv(cache::BasicILMCache) = _y_griddiv(cache,Val(cache_datatype(cache)))

_y_griddiv(cache,::Val{:scalar}) = nothing

function _y_griddiv(cache,::Val{:vector})
    _,yc = coordinates(cache.gdiv_cache,cache.g)
    p = zeros_griddiv(cache)
    p .= yc'
end

"""
    y_gridgrad(::BasicILMCache)

Return basic grid gradient field data filled with the grid `y` coordinate.
If the grid data is scalar, then the output of this function is
of `Edges{Primal}` type and the coordinate field for each component
can be accessed with the `u` and `v` field, respectively. If the grid
data is vector, then the output is of type `EdgeGradient{Primal}`,
and the coordinates are in `dudx`, `dudy`, `dvdx`, and `dvdy` fields.
"""
function y_gridgrad(cache::BasicILMCache{N,SCA,GT}) where {N,SCA,GT<:Nodes{Primal}}
    p = zeros_gridgrad(cache)
    _, yu = coordinates(p.u,cache.g)
    _, yv = coordinates(p.v,cache.g)
    p.u .= yu'
    p.v .= yv'
    return p
end

function y_gridgrad(cache::BasicILMCache{N,SCA,GT}) where {N,SCA,GT<:Edges{Primal}}
    p = zeros_gridgrad(cache)
    _, ydudx = coordinates(p.dudx,cache.g)
    _, ydudy = coordinates(p.dudy,cache.g)
    p.dudx .= ydudx'
    p.dudy .= ydudy'
    p.dvdx .= ydudy'
    p.dvdy .= ydudx'
    return p
end

# Extend operators on body points
"""
    points(cache::BasicILMCache)

Return the coordinates (as `VectorData`) of the surface points associated with `cache`
"""
points(cache::BasicILMCache) = points(cache.bl)

"""
    points(cache::PointCollectionCache)

Return the coordinates (as `VectorData`) of the surface points associated with `cache`
"""
points(cache::PointCollectionCache) = cache.pts

"""
    normals(cache::BasicILMCache)

Return the normals (as `VectorData`) of the surface points associated with `cache`
"""
normals(cache::BasicILMCache) = cache.nrm

"""
    areas(cache::BasicILMCache)

Return the areas (as `ScalarData`) of the surface panels associated with `cache`
"""
areas(cache::BasicILMCache) = cache.ds

# Extend norms and inner products
"""
    norm(u::GridData,cache::BasicILMCache)

Calculate the norm of grid data `u`, using the scaling associated with `cache`.
""" norm(u::GridData,cache::BasicILMCache)

norm(u::GridData,cache::BasicILMCache{N,GridScaling}) where {N} = norm(u,cache.g)

norm(u::GridData,cache::BasicILMCache{N,IndexScaling}) where {N} = norm(u)

"""
    dot(u1::GridData,u2::GridData,cache::BasicILMCache)

Calculate the inner product of grid data `u1` and `u2`, using the scaling associated with `cache`.
""" dot(u1::GridData,u2::GridData,cache::BasicILMCache)

dot(u1::GridData,u2::GridData,cache::BasicILMCache{N,GridScaling}) where {N} = dot(u1,u2,cache.g)

dot(u1::GridData,u2::GridData,cache::BasicILMCache{N,IndexScaling}) where {N} = dot(u1,u2)

"""
    norm(u::PointData,cache::BasicILMCache)

Calculate the norm of surface point data `u`, using the scaling associated with `cache`.
""" norm(u::PointData,cache::BasicILMCache)

norm(u::PointData{N},cache::BasicILMCache{N,GridScaling}) where {N} = norm(u,cache.ds)

norm(u::PointData{N},cache::BasicILMCache{N,IndexScaling}) where {N} = norm(u)

"""
    norm(u::PointData,cache::BasicILMCache,i::Int)

Calculate the norm of surface point data `u`, using the scaling associated with `cache`,
for body `i` in the body list of `cache`.
""" norm(u::PointData,cache::BasicILMCache,i::Int)

norm(u::PointData{N},cache::BasicILMCache{N,GridScaling},i::Int) where {N} = norm(u,cache.ds,cache.bl,i)

norm(u::PointData{N},cache::BasicILMCache{N,IndexScaling},i::Int) where {N} = norm(u,cache.bl,i)

for f in [:ScalarData,:VectorData]
  @eval dot(u1::$f{N},u2::$f{N},cache::BasicILMCache{N,GridScaling}) where {N} = dot(u1,u2,cache.ds)

  @eval dot(u1::$f{N},u2::$f{N},cache::BasicILMCache{N,IndexScaling}) where {N} = dot(u1,u2)
end
"""
    dot(u1::PointData,u2::PointData,cache::BasicILMCache)

Calculate the inner product of surface point data `u1` and `u2`, using the scaling associated with `cache`.
""" dot(u1::PointData,u2::PointData,cache::BasicILMCache)

"""
    dot(u1::PointData,u2::PointData,cache::BasicILMCache,i)

Calculate the inner product of surface point data `u1` and `u2` for body
`i` in the cache `cache`, scaling as appropriate for this cache.
""" dot(u1::PointData,u2::PointData,cache::BasicILMCache,i::Int)

dot(u1::PointData{N},u2::PointData{N},cache::BasicILMCache{N,GridScaling},i::Int) where {N} = dot(u1,u2,cache.ds,cache.bl,i)

dot(u1::PointData{N},u2::PointData{N},cache::BasicILMCache{N,IndexScaling},i::Int) where {N} = dot(u1,u2,cache.bl,i)

## Integration
"""
    integrate(u::PointData,cache::BasicILMCache)

Calculate the discrete surface integral of data `u` on the immersed points in `cache`.
This uses trapezoidal rule quadrature. If `u` is `VectorData`, then this returns a vector of the integrals in
each coordinate direction. This operation produces the same effect,
regardless if `cache` is set up for `GridScaling` or `IndexScaling`. In both
cases, the surface element areas are used.
"""
@inline integrate(u::PointData{N},cache::BasicILMCache{N}) where {N} = integrate(u,cache.ds)

"""
    integrate(u::PointData,cache::BasicILMCache,i::Int)

Calculate the discrete surface integral of scalar data `u` on the immersed points in `cache`,
on body `i` in the body list in `cache`.
This uses trapezoidal rule quadrature. If `u` is `VectorData`, then this returns a vector of the integrals in
each coordinate direction. This operation produces the same effect,
regardless if `cache` is set up for `GridScaling` or `IndexScaling`. In both
cases, the surface element areas are used.
"""
@inline integrate(u::PointData{N},cache::BasicILMCache{N},i::Int) where {N} = integrate(u,cache.ds,cache.bl,i)

# Extending some operations on body lists to the enclosing cache
"""
    view(u::PointData,cache::BasicILMCache,i::Int)

Provide a `view` of point data `u` corresponding to body `i` in the
list of bodies in `cache`.
"""
view(u::PointData,cache::BasicILMCache,i::Int) = view(u,cache.bl,i)

"""
    copyto!(u::PointData,v::PointData,cache::BasicILMCache,i::Int)

Copy the data in the elements of `v` associated with body `i` in the body list in `cache` to
the corresponding elements in `u`. These data must be of the same type (e.g.,
`ScalarData` or `VectorData`) and have the same length.
""" copyto!(u::PointData,v::PointData,cache::BasicILMCache,i::Int)

"""
    copyto!(u::ScalarData,v::AbstractVector,cache::BasicILMCache,i::Int)

Copy the data in `v` to the elements in `u` associated with body `i` in the body list in `cache`.
`v` must have the same length as this subarray of `u` associated with `i`.
""" copyto!(u::ScalarData,v::AbstractVector,cache::BasicILMCache,i::Int)

@inline copyto!(u::PointData{N},v,cache::BasicILMCache{N},i::Int) where {N} = copyto!(u,v,cache.bl,i)

"""
    arcs(cache::BasicILMCache)

Return `ScalarData` of arc length coordinates for the body surface(s) in
the given cache `cache`.
"""
arcs(cache::BasicILMCache) = arcs(cache.bl)


text:
using Base.Meta

"""
Finds ISPC kernels in `func`, extracts them, and
replaces them with calls to the compiled ISPC functions.
"""
function make_ispc_main(func::Expr)
    substitute(func) do expr
        if isa(expr, Expr) && expr.head == :call
            func = expr.args[1]
            if func == GlobalRef(ISPC, :ispc_kernel)
                identifier, lambda, options = expr.args[2:end]

                ast = get_ast(lambda)
                arg_names = ast.args[1]
                captured = ast.args[2][2]
                body = ast.args[3]

                # Find out which arguments are modified by the
                # kernel call:
                modified = []
                kernel_argnames = Any[arg_names...]
                for (name, typ, flags) in captured
                    push!(kernel_argnames, name)
                    if (flags & 0x4) != 0
                        # Fortunately, this is already set after lowering,
                        # before type inference.
                        push!(modified, name)
                    end
                end

                println("Extracted kernel $identifier$((kernel_argnames...)) -- modified: $((modified...))")
                # println(strip_lineno(body))

                # identifier is a Val{symbol}, unbox the symbol:
                id = identifier.parameters[1]

                # Register the fragment:
                ISPC.ispc_fragments[id] = lambda
                ISPC.ispc_fragment_opts[id] = options

                # Generate the corresponding call.
                 # Don't forget we're dealing with lowered ASTs, expand:
                kernel_call = expand(:(ISPC.kernel_call($identifier,
                                             $((kernel_argnames...,)))))

                # Don't forget to handle any returned values:
                if isempty(modified)
                    return kernel_call
                else
                    return :(
                        $(modified...) = $kernel_call
                    )
                end
            end
        end
        nothing # not matched, don't substitute
    end
end

"""
Finds ranges of statements inside opening and closing ISPC tags.
"""
function ispc_ranges(statements)
    ranges = []
    cur_id = nothing
    start = 1
    for (i, stmt) in enumerate(statements)
        if isa(stmt, Expr) && stmt.head == :meta
            if stmt.args[1] == :ispc
                identifier = stmt.args[2]
                if cur_id == nothing && identifier != :coherent
                    # opening tag
                    prev_range = start:i-1
                    if !isempty(prev_range)
                        push!(ranges, (prev_range, nothing))
                    end
                    cur_id = identifier
                    start = i
                elseif cur_id == identifier
                    # closing tag
                    cur_id = nothing
                    push!(ranges, (start:i, identifier))
                    start = i+1
                end
            end
        end
    end
    @assert cur_id == nothing
    tail_range = start:length(statements)
    if !isempty(tail_range)
        push!(ranges, (tail_range, nothing))
    end
    return ranges
end


text:
"""
Takes `frames`, which is supposed to be an array of images,
saves them as png's at path and then creates an webm video
from that with the name `name`
"""
function create_video(frames::Vector, name, screencap_folder, resample_steps=0, remove_destination=true)
    println("saving frames for $name")
    mktempdir() do path
        frame1 = first(frames)
        for i=1:resample_steps
            frame1 = Images.restrict(frame1)
        end
        resolution = size(frame1)
        for (i,frame) in enumerate(frames)
            resampled = frame
            for x=1:resample_steps
                resampled = Images.restrict(resampled)
            end
            frame = map(RGB{U8}, resampled)
            save(joinpath(path, "$name$i.png"), frame, true)
        end
        len = length(frames)
        frames = [] # free frames...
        oldpath = pwd()
        cd(path)
        mktemp(path) do io, path
            # output stdout to tmpfile , since there is too much going on with it
            run(pipeline(
                `png2yuv -I p -f 30 -b 1 -n $len -j $name%d.png`,
                stdout="$(name).yuv",stdin=io, stderr=io
            ))
            run(pipeline(
                `vpxenc --good --cpu-used=0
                    --auto-alt-ref=1 --lag-in-frames=16 --end-usage=vbr
                    --passes=1 --threads=4 --target-bitrate=3500
                    -o $(name).webm $(name).yuv`,
                stdout=io, stdin=io, stderr=io
            ))
        end
        targetpath = abspath(joinpath(screencap_folder, "$(name).webm"))
        sourcepath = abspath("$(name).webm")
        mv(sourcepath, targetpath, remove_destination=remove_destination)
        cd(oldpath)
    end
end

function create_video(frame, name, path, resample_steps=0)
    println("saving static image for $name")
    targetpath = abspath(joinpath(path, "$(name).png"))
    for i=1:resample_steps
        frame = Images.restrict(frame)
    end
    save(targetpath, frame, true)
end
#downsample
#`ffmpeg -i volume.webm -codec:v libvpx -quality good -cpu-used 0 -b:v 500k -qmin 10 -qmax 42 -maxrate 500k -bufsize 1000k -threads 4 -vf scale=-1:480 -codec:a libvorbis -b:a 128k volume2.webm`


text:
precompile(SphericalHarmonicsCache, (Type{Float64}, Int, Type{FullRange}, ComplexHarmonics))
precompile(SphericalHarmonicsCache, (Type{Float64}, Int, Type{FullRange}, RealHarmonics))
precompile(SphericalHarmonicsCache, (Type{Float64}, Int, Type{ZeroTo}, ComplexHarmonics))
precompile(SphericalHarmonicsCache, (Type{Float64}, Int, Type{ZeroTo}, RealHarmonics))
_S = SphericalHarmonicsCache(Float64, 0, FullRange, ComplexHarmonics())
precompile(computePlmx!, (typeof(_S), Float64, Int))
precompile(computePlmx!, (typeof(_S.P), Float64, Int))
precompile(computePlmcostheta!, (typeof(_S), Float64, Int))
precompile(computePlmcostheta!, (typeof(_S.P), Float64, Int))
precompile(computeYlm!, (typeof(_S), Float64, Float64, Int))


text:
push!(LOAD_PATH,"../src/")

using Documenter
using ScikitSpatial

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

makedocs(
    sitename="ScikitSpatial.jl",
    modules=[ScikitSpatial],
    doctest=true,
    pages = [
        "Introduction" => "index.md",
        "Types" => [
            "types/line.md",
            "types/plane.md",
        ],
        "Computations" => [
            "Measurement" => "computations/measurement.md",
            "Comparison" => "computations/comparison.md",
            "Projection" => "computations/projection.md",
            "Intersection" => "computations/intersection.md",
            "Transformation" => "computations/transformation.md",
        ],
    ]
)

deploydocs(
    repo = "github.com/ajhynes7/ScikitSpatial.jl.git",
)


text:
# Code from @devmotion
# https://github.com/devmotion/\
# CalibrationErrorsDistributions.jl/blob/main/src/distances/bures.jl
using ExactOptimalTransport

using LinearAlgebra
using Random
using PDMats

@testset "bures.jl" begin
    function _sqbures(A, B)
        sqrt_A = sqrt(A)
        return tr(A) + tr(B) - 2 * tr(sqrt(sqrt_A * B * sqrt_A'))
    end

    function rand_matrices(n)
        A = randn(n, n)
        B = A' * A + I
        return B, PDMat(B), PDiagMat(diag(B)), ScalMat(n, B[1])
    end

    for (x, y) in Iterators.product(rand_matrices(10), rand_matrices(10))
        xfull = Matrix(x)
        yfull = Matrix(y)
        @test ExactOptimalTransport.sqbures(x, y) ≈ _sqbures(xfull, yfull)
    end
end


text:
using ADCME
using PyCall
using LinearAlgebra
using PyPlot
using Random
Random.seed!(233)

function mat_vec(a,x)
    mat_vec_ = load_op_and_grad("./build/libMatVec","mat_vec")
    a,x = convert_to_tensor(Any[a,x], [Float64,Float64])
    mat_vec_(a,x)
end

# TODO: specify your input parameters
a0 = rand(10,10)
x0 = rand(10)
a = constant(a0)
u = mat_vec(a0,x0)
sess = Session(); init(sess)
@show run(sess, u)≈a0 * x0