text:
"""
    NullCorrection
Implements Castaing's First Approximation (i.e. No correction Z⋅A = 1)
"""
struct NullCorrection <: MatrixCorrection
    material::Material
    subshell::AtomicSubShell
    E0::AbstractFloat

    NullCorrection(mat::Material, ashell::AtomicSubShell, e0) = new(mat, ashell, e0)
end

Base.show(io::IO, nc::NullCorrection) = print(io, "Null[ $(nc.material), $(nc.subshell), $(0.001 * nc.E0) keV]")

F(mc::NullCorrection) = 1.0
Fχ(mc::NullCorrection, xray::CharXRay, θtoa::Real) = 1.0
NeXLCore.atomicsubshell(mc::NullCorrection) = mc.subshell
NeXLCore.material(mc::NullCorrection) = mc.material
beamEnergy(mc::NullCorrection) = mc.E0 # in eV
ϕ(mc::NullCorrection, ρz) = @assert "Does not make sense for the null correction."
ϕabs(mc::MatrixCorrection, ρz, θtoa) = @assert "Does not make sense for the null correction."
NeXLMatrixCorrection.matrixcorrection(::Type{NullCorrection}, mat::Material, ass::AtomicSubShell, e0::Float64) = NullCorrection(mat,ass,e0)

max_stars_repo_name: NicholasWMRitchie/NeXLMatrixCorrection

max_issues_repo_name: NicholasWMRitchie/NeXLMatrixCorrection

max_forks_repo_name: NicholasWMRitchie/NeXLMatrixCorrection


text:
using FileIO, BedgraphFiles
using Bedgraph

using IteratorInterfaceExtensions
using TableTraits

using DataFrames
using Query

using Test
using Logging

# old_logger = global_logger(ConsoleLogger(stdout, Logging.Debug))

module Bag

using Bedgraph

const chroms = ["chr19", "chr19", "chr19", "chr19", "chr19", "chr19", "chr19", "chr19", "chr19"]
const firsts = [49302000, 49302300, 49302600, 49302900, 49303200, 49303500, 49303800, 49304100, 49304400]
const lasts = [49302300, 49302600, 49302900, 49303200, 49303500, 49303800, 49304100, 49304400, 49304700]
const values = [-1.0, -0.75, -0.50, -0.25, 0.0, 0.25, 0.50, 0.75, 1.00]

const record = Bedgraph.Record("chr1", 1, 1, 0)
const records = Bedgraph.Record.(Bag.chroms, Bag.firsts, Bag.lasts, Bag.values)

const file = joinpath(@__DIR__, "data.bedgraph")
const file_headerless = joinpath(@__DIR__, "data-headerless.bedgraph")

const tmp_output_path = tempname() * ".bedgraph"

end # module Bag

using .Bag

@testset "BedgraphFiles" begin

@test isfile(Bag.file) == true
@test isfile(Bag.file_headerless) == true

# Load tests.
loader = load(Bag.file)
@test IteratorInterfaceExtensions.isiterable(loader) == true
@test TableTraits.isiterabletable(loader) == true

loaded = Vector{Bedgraph.Record}(loader)
@test Vector{Bedgraph.Record} == typeof(loaded)

loader_from_headerless = load(Bag.file_headerless)
@test isiterable(loader_from_headerless) == true
@test TableTraits.isiterabletable(loader_from_headerless) == true

loaded_from_headerless = Vector{Bedgraph.Record}(loader_from_headerless)
@test Vector{Bedgraph.Record} == typeof(loaded_from_headerless)

@test IteratorInterfaceExtensions.isiterable(Bag.records) == true
@test TableTraits.isiterabletable(Bag.records) == true

@test Bag.records == loaded
@test Bag.records == loaded_from_headerless

# Save and load from Vector{Bedgraph.Record}.
save(Bag.tmp_output_path, Bag.records)

@debug "direct load into Vector{Bedgraph.Record} - commencing"
@test Bag.records == Vector{Bedgraph.Record}(load(Bag.tmp_output_path))
@debug "direct load into Vector{Bedgraph.Record} - complete"

@test Bag.records == load(Bag.tmp_output_path) |> Vector{Bedgraph.Record}

# Save usign query.
Bag.records |> save(Bag.tmp_output_path)
@test Bag.records == Vector{Bedgraph.Record}(load(Bag.tmp_output_path))
@test Bag.records == load(Bag.tmp_output_path) |> Vector{Bedgraph.Record}

# Check return of data from save method.
@test Bag.records == Bag.records |> save(Bag.tmp_output_path)

# Check piping/continuations through Query.jl.
load("data.bedgraph") |> @filter(_.chrom == "chr19" && _.first > 49302900 && _.last < 49303800) |> save(Bag.tmp_output_path)
@test [Bedgraph.Record("chr19", 49303200, 49303500, 0.0)] == load(Bag.tmp_output_path) |> Vector{Bedgraph.Record}

@testset "Integrations" begin
    include("integrations/test-DataFrames.jl")
end # testset Transformers

println()
show(load(Bag.file))
println()

end

ext: jl


text:
using ADI
using CoordinateTransformations
using Statistics
using WoodburyMatrices
using LinearAlgebra
using HCIToolbox
using Interpolations
using PSFModels: PSFModel
using StaticArrays

using ImageTransformations: box_extrapolation, center
using Interpolations: AbstractExtrapolation

import Distributions
import Distributions: ContinuousMatrixDistribution

struct JointModel{T <: AbstractArray, M <: AbstractMatrix, V <: AbstractVector, W <: AbstractWoodbury,WT,G<:CubeGenerator}
    cube::T
    target::M
    angles::V
    Σ::W
    Σ⁻¹AΛAᵀ::M
    w̄A::WT
    gen::G
    tmpR::M
end

"""
    JointModel(A, w, cube, psf)

Construct a joint model using the given basis vectors `A` and weights `w`. They should have shape `(n, Npx)` and `(n, M)`, respectively. `cube` will be used as the target for the likelihood.

# Extended help

The mathematical formulation of this model is
```math
\\text{cube} \\sim \\mathbf{A} \\cdot \\mathbf{w} + \\text{PSF} + \\epsilon
```
"""
function JointModel(A, w, cube::AbstractArray{T,3}, angles, psf; kwargs...) where T
    w̄ = mean(w, dims=1)
    w̄A = w̄ * A
    X = flatten(cube)
    resid = X - w * A
    Λ = cov(w, dims=1)
    C = fake_covariance(resid)
    Σ = SymWoodbury(C, A', Λ)
    Σ⁻¹AΛAᵀ = Σ \ (A' * (Λ * A))
    target = X .- w̄A
    gen = CubeGenerator(cube, angles, psf; kwargs...)
    return JointModel(cube, target, angles, Σ, Σ⁻¹AΛAᵀ, w̄A, gen, similar(target))
end

function JointModel(A, w, cube::AnnulusView{T,3}, angles, psf; kwargs...) where T
    w̄ = mean(w, dims=1)
    w̄A = w̄ * A
    X = cube()
    resid = X - w * A
    Λ = cov(w, dims=1)
    C = fake_covariance(resid)
    Σ = SymWoodbury(C, A', Λ)
    Σ⁻¹AΛAᵀ = Σ \ (A' * (Λ * A))
    target = X .- w̄A
    gen = CubeGenerator(cube, angles, psf; kwargs...)
    return JointModel(cube, target, angles, Σ, Σ⁻¹AΛAᵀ, w̄A, gen, similar(target))
end

"""
    Distributions.loglikelihood(::JointModel, [data]; r, theta, A=1)
    Distributions.loglikelihood(::JointModel, [data]; x, y, A=1)

Return the likelihood of the model target given the input parameters. By default uses the model's cube as the input, but a matrix can be passed directly. Returns the logkernel of a matrix normal distribution without column-wise variance.
"""
function Distributions.loglikelihood(model::JointModel, pos; A=zero(eltype(model.cube)))
    T = float(typeof(A))
    base = T.(model.target)
    R = model.gen(base, pos; A=-A)
    return -tr(R * (model.Σ \ R')) / 2
end

function ADI.reconstruct(model::JointModel, pos; A=zero(eltype(model.cube)))
    T = float(typeof(A))
    base = T.(model.target)
    R = model.gen(base, pos; A=-A)
    return expand(model.w̄A .+ R * model.Σ⁻¹AΛAᵀ)
end

function ADI.reconstruct(model::JointModel{<:AnnulusView}, pos; A=zero(eltype(model.cube)))
    T = float(typeof(A))
    base = T.(model.target)
    R = model.gen(base, pos; A=-A)
    return inverse(model.cube, model.w̄A .+ R * model.Σ⁻¹AΛAᵀ)
end

function model(model::JointModel, pos; A=zero(eltype(model.cube)))
    T = float(typeof(A))
    return model.gen(T, pos; A)
end

function Statistics.mean(m::JointModel, args...; params...)
    plan = model(m, args...; params...)
    L = reconstruct(m, args...; params...)
    return expand(plan .+ L)
end

function Statistics.mean(m::JointModel{<:AnnulusView}, args...; params...)
    plan = model(m, args...; params...)
    L = reconstruct(m, args...; params...)
    return inverse(m.cube, plan .+ L)
end

# TODO
# probably use Distributions.logkernel and Distributions.loglikelihood to further
# allow this to become a "distribution" for use in other inference methods.
struct JointDistribution{T,M<:AbstractMatrix{T},W<:AbstractWoodbury} <: ContinuousMatrixDistribution
    μ::M
    Σ::W
end

function Distributions.loglikelihood(d::JointDistribution, X::AbstractMatrix)
    R = X - d.μ
    return -tr(R * (d.Σ \ R')) / 2
end

Base.size(d::JointDistribution) = size(d.μ)

# covariance is just a diagonal of the variance
function fake_covariance(mat::AbstractMatrix)
    D = var(mat, dims = 1) |> vec
    @. D[iszero(D)] = Inf
    return Diagonal(D)
end


text:
using DashBootstrapComponents, DashHtmlComponents

offcanvas = html_div([
    dbc_button("Open scrollable offcanvas", id = "open-offcanvas-scrollable", n_clicks = 0),
    dbc_offcanvas(
        html_p("The contents on the main page are now scrollable."),
        id = "offcanvas-scrollable",
        scrollable = true,
        title = "Scrollable Offcanvas",
        is_open = false,
    ),
])

callback!(
    app,
    Output("offcanvas-scrollable", "is_open"),
    Input("open-offcanvas-scrollable", "n_clicks"),
    State("offcanvas-scrollable", "is_open"),
) do n1, is_open
    return n1 > 0 ? is_open == 0 : is_open
end;

max_stars_repo_name: glsdown/dash-bootstrap-components

max_issues_repo_name: glsdown/dash-bootstrap-components

max_forks_repo_name: glsdown/dash-bootstrap-components


text:
# Toy problem to test the solver

using luteos

P = P3() # Polynomial order of solution

mesh   = Mesh2D( "square", P, N = 30, M = 30)
master = Master2D( P )

# source function
function funcS( p::Array{Float64} )
  return fill(1.0, size(p,1), 1)
end

function funcB( p::Array{Float64} )
  return fill(0.0, size(p,1), 1)
end

prob = CDR( "Example 3 - Poission", 1.0, [0.0,0.0], funcS, [1,1,1,1], true, [funcB, funcB, funcB, funcB] )

@time (uhath, uh, qh, uhathTri ) = hdgSolve( master, mesh, prob )

# write solution
writeTecplot( "blaCD.dat", prob, mesh, uh, qh )


text:
module SNNTest
using LightGraphs
using SpikingNeuralNets, SpikingNeuralNets.LIFs, SpikingNeuralNets.SNNPlot

const iterations = 1000
const dt = 1.0 # timestep

I(snn, n, s) = 1.00
s(snn, n, s) = s

graph = Graph(1)
add_edge!(graph, 1, 1)
lif = LIF{Float64, Int}(graph; V=[[0.0]], dt=dt, Vrest=[0.0], Vθ=[1.5],
                        Vreset=[0.0], τref=[10.0], τ=[100.0], R=[1.75])
add_channel!(lif, "input", I, s, [Float64[]])
lif, V = step!(lif, iterations, voltages; threaded=false)
V = reduce(hcat, V) # convert V from a Vector{Vector} to a Matrix
S = spiketrain(lif, iterations)

# Make a raster plot of the output
display(rasterplot(S; xlim=[0, iterations], timestep=dt))
display(voltageplot(V; xlim=[0, iterations], timestep=dt))

end


text:
@testset "IntVar.jl" begin
    @testset "isbound()" begin
        trailer = SeaPearl.Trailer()
        x = SeaPearl.IntVar(2, 6, "x", trailer)

        @test !SeaPearl.isbound(x)

        SeaPearl.assign!(x, 3)

        @test SeaPearl.isbound(x)
    end

    @testset "assignedValue()" begin
        trailer = SeaPearl.Trailer()
        x = SeaPearl.IntVar(5, 5, "x", trailer)

        @test SeaPearl.assignedValue(x) == 5

        y = SeaPearl.IntVar(5, 8, "y", trailer)
        @test_throws AssertionError SeaPearl.assignedValue(y)
    end

    @testset "IntVar()" begin
        trailer = SeaPearl.Trailer()
        x = SeaPearl.IntVar(5, 8, "x", trailer)

        @test length(x.domain) == 4
        @test 5 in x.domain && 8 in x.domain && !(4 in x.domain) && !(9 in x.domain)
    end

    @testset "assign!()" begin
        trailer = SeaPearl.Trailer()
        x = SeaPearl.IntVar(5, 6, "x", trailer)

        @test_throws AssertionError SeaPearl.assignedValue(x)

        SeaPearl.assign!(x, 5)

        @test SeaPearl.assignedValue(x) == 5
    end
end


text:
#===============================================================================

Using a combination of black square tiles and oblong tiles chosen from: red
tiles measuring two units, green tiles measuring three units, and blue tiles
measuring four units, it is possible to tile a row measuring five units in
length in exactly fifteen different ways.

How many ways can a row measuring fifty units in length be tiled?

===============================================================================#

using Memoize
using Base.Test

@memoize function ncases_swb(rowlen)
    # Return the number of cases starting with a block

    ncases_swglb(2, rowlen) + ncases_swglb(3, rowlen) + ncases_swglb(4, rowlen)
end

@memoize function ncases_swglb(blocklen, rowlen)
    # Return the number of cases starting with a block of a given length
    if rowlen < blocklen
        return 0
    end

    result = 1
    # case for which there is only the starting block

    rowlen -= blocklen
    # account for the head block

    for nspace = 0:rowlen-2
        # nbr spaces allowed between first block and second one
        result += ncases_swb(rowlen-nspace)
    end

    result
end

function ncases(rowlen)
    result = 1
    # No block at all

    for nspace = 0:rowlen-2
        # nbr spaces before first block
        result += ncases_swb(rowlen-nspace)
    end

    result
end

@test ncases(5) == 15

@printf("Euler 0117: %d", ncases(50))

max_stars_repo_name: dpieroux/euler

max_issues_repo_name: dpieroux/euler

max_forks_repo_name: dpieroux/euler


text:
# ---
# title: 1646. Get Maximum in Generated Array
# id: problem1646
# author: Tian Jun
# date: 2020-10-31
# difficulty: Easy
# categories: Array
# link: <https://leetcode.com/problems/get-maximum-in-generated-array/description/>
# hidden: true
# ---
#
# You are given an integer `n`. An array `nums` of length `n + 1` is generated
# in the following way:
#
#   * `nums[0] = 0`
#   * `nums[1] = 1`
#   * `nums[2 * i] = nums[i]` when `2 <= 2 * i <= n`
#   * `nums[2 * i + 1] = nums[i] + nums[i + 1]` when `2 <= 2 * i + 1 <= n`
#
# Return ****_the **maximum** integer in the array _`nums`.
#
#
#
# **Example 1:**
#
#
#
#     Input: n = 7
#     Output: 3
#     Explanation: According to the given rules:
#       nums[0] = 0
#       nums[1] = 1
#       nums[(1 * 2) = 2] = nums[1] = 1
#       nums[(1 * 2) + 1 = 3] = nums[1] + nums[2] = 1 + 1 = 2
#       nums[(2 * 2) = 4] = nums[2] = 1
#       nums[(2 * 2) + 1 = 5] = nums[2] + nums[3] = 1 + 2 = 3
#       nums[(3 * 2) = 6] = nums[3] = 2
#       nums[(3 * 2) + 1 = 7] = nums[3] + nums[4] = 2 + 1 = 3
#     Hence, nums = [0,1,1,2,1,3,2,3], and the maximum is 3.
#
#
# **Example 2:**
#
#
#
#     Input: n = 2
#     Output: 1
#     Explanation: According to the given rules, the maximum between nums[0], nums[1], and nums[2] is 1.
#
#
# **Example 3:**
#
#
#
#     Input: n = 3
#     Output: 2
#     Explanation: According to the given rules, the maximum between nums[0], nums[1], nums[2], and nums[3] is 2.
#
#
#
#
# **Constraints:**
#
#   * `0 <= n <= 100`
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end


text:
try
    using RCall
catch
end
using RCall
using Distributions
using StatsBase
using Iterators
using Gurobi
using JLD

using EmpiricalBayes

nreps =20
k_max = 12 #in total 25*16 replications

true_dists = [MixtureModel([ Normal(-0.3,.2), Normal(0,0.9)],[0.8, 0.2]), EmpiricalBayes.ash_flattop ]
xs = collect(-3:0.5:3)
minimax_σs = [0.2, 0.5]
ms = [20_000, 100_000]

combs = []

for (true_dist, x, σ, m) in product(true_dists, xs, minimax_σs, ms)
    push!(combs, Dict(:dist=>true_dist, :x=>x, :σ=>σ, :m=>m) )
end

i = parse(Int64, ARGS[1])
i=1
comb = combs[i]

m = comb[:m]
σ = comb[:σ]
true_dist = comb[:dist]
x = comb[:x]

marginal_grid = collect(range(-7.0, stop=7.0, length=1001));
prior_grid = collect(range(-3, stop=3, length=121));
marginal_h = marginal_grid[2]-marginal_grid[1];

d_true = NormalConvolutionProblem(true_dist, marginal_grid)
ds = MixingNormalConvolutionProblem(Normal, σ, prior_grid, marginal_grid);

target = PosteriorTarget(LFSRNumerator(x))

dry_run = CEB_ci(rand(d_true,Int(m/2)), rand(d_true, Int(m/2)), ds, target)

dry_run_snoop = CEB_ci(rand(d_true,m), rand(d_true, m), ds, target)

for k=1:k_max
    sim_array = []
    for l = 1:nreps
        @show l
        Xs = rand(d_true,m)
        brad_non_null =  fit(BradDeconvolveR, Xs; deltaAt = :nothing)
        brad_null =  fit(BradDeconvolveR, Xs)

        n_total = length(Xs)
        n_half = ceil(Int, n_total/2)
        idx_test = sample(1:n_total, n_half, replace=false)
        idx_train = setdiff(1:n_total, idx_test)
        Xs_train = Xs[idx_train]
        Xs_test = Xs[idx_test]

        ceb_res = CEB_ci(Xs_train, Xs_test, ds, target,
                    dry_run[2], dry_run[3], dry_run[4].C_bias)
        ceb_res_snoop = CEB_ci(Xs, Xs, ds, target,
                    dry_run_snoop[2], dry_run_snoop[3], dry_run_snoop[4].C_bias)

        push!(sim_array, (ceb_res, ceb_res_snoop, brad_null, brad_non_null))
    end
    res = (i, comb, sim_array)
    save("/scratch/users/ignat/sims/May22/mysim_$(i)_$(k).jld", "res", res)
end


text:
"""
    CubicEoSDatabase.NotFoundError(what, where) <: Exception

Describes `where` `what` was not found. The `where` should be database object.
"""
struct NotFoundError <: Exception
    what
    where
end

Base.show(io::IO, x::NotFoundError) = begin
    print(typeof(x), "\nKeyword $(repr(x.what)) not found in $(x.where)")
end


text:
module KrigingModel

# using M:f  = call `f` only
# import M   = call `M.anyfunction`
# import M:f = extend `f`

using Statistics: mean!, mean, std
using ElasticArrays: ElasticArray

import GaussianProcesses, LinearAlgebra

const GP = GaussianProcesses

include("scaler/scaler.jl")
include("optimizer/optimizer.jl")
include("Kriging/Kriging.jl")

export StandardScaler, MinMaxScaler
export CRSOptimizer, ISRESOptimizer, LBFGSOptimizer, BOBYQAOptimizer
export minimize, maximize
export Kriging, tune!, predict_full, getx, gety, get_samples

end # module


text:
module ReverseDiffExpressions

# using VectorizationBase, SIMDPirates
# using StackPointers, LinearAlgebra, PaddedMatrices#, StructuredMatrices, ProbabilityDistributions, Parameters
using LoopVectorization, UnPack, OffsetArrays
using LoopVectorization: Instruction, LoopSet, instruction
using ReverseDiffExpressionsBase: Target, InstructionArgs, DiffRule, DERIVATIVERULES, FALLBACK_RULES, num_parents
using ReverseDiffExpressionsBase

export Target

# using PaddedMatrices:
    # AbstractFixedSizeVector, AbstractMutableFixedSizeVector, UninitializedVector,
    # AbstractFixedSizeMatrix, AbstractMutableFixedSizeMatrix, UninitializedMatrix,
    # AbstractFixedSizeArray, AbstractMutableFixedSizeArray, UninitializedArray

# using StructuredMatrices: AbstractLowerTriangularMatrix
# using ProbabilityDistributions: distribution_diff_rule!

# using FunctionWrappers: FunctionWrapper
# using MacroTools: postwalk, prewalk, @capture, @q

# import SIMDPirates: vsum, vadd, vifelse
# import ReverseDiffExpressionsBase: adj, ∂mul, ∂getindex, RESERVED_INCREMENT_SEED_RESERVED!,
    # InitializedVarTracker, add_aliases!, allocate!, initialize!, isallocated, deallocate!

# FIXME: create it; stub for now

const STACK_POINTER_NAME = Symbol("###STACK##POINTER###")

include("variables.jl")
include("funcs.jl")
include("models.jl")
include("tracking_description.jl")
include("differentiate.jl")
include("loopsets/loopset_gradients.jl") # defines a module that includes the other loopsets/*.jl files
include("diffrulevars.jl")
include("lowering.jl")
# include("adjoints.jl")
# include("misc_functions.jl")
# include("special_diff_rules.jl")
# include("reverse_autodiff_passes.jl")

# include("precompile.jl")
# _precompile_()

# @def_stackpointer_fallback emax_dose_response ITPExpectedValue ∂ITPExpectedValue HierarchicalCentering ∂HierarchicalCentering
function __init__()
    # @add_stackpointer_method emax_dose_response ITPExpectedValue ∂ITPExpectedValue HierarchicalCentering ∂HierarchicalCentering
end

end # module


text:
# Autogenerated wrapper script for Lua_jll for i686-w64-mingw32
export liblua, lua, luac

using Readline_jll
JLLWrappers.@generate_wrapper_header("Lua")
JLLWrappers.@declare_library_product(liblua, "lua54.dll")
JLLWrappers.@declare_executable_product(lua)
JLLWrappers.@declare_executable_product(luac)
function __init__()
    JLLWrappers.@generate_init_header(Readline_jll)
    JLLWrappers.@init_library_product(
        liblua,
        "bin\\lua54.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_executable_product(
        lua,
        "bin\\lua.exe",
    )

    JLLWrappers.@init_executable_product(
        luac,
        "bin\\luac.exe",
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


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

mutable struct VirtualMachineScaleSetVMExtensionsSummary <: SwaggerModel
    name::Any # spec type: Union{ Nothing, String } # spec name: name
    statusesSummary::Any # spec type: Union{ Nothing, Vector{VirtualMachineStatusCodeCount} } # spec name: statusesSummary

    function VirtualMachineScaleSetVMExtensionsSummary(;name=nothing, statusesSummary=nothing)
        o = new()
        validate_property(VirtualMachineScaleSetVMExtensionsSummary, Symbol("name"), name)
        setfield!(o, Symbol("name"), name)
        validate_property(VirtualMachineScaleSetVMExtensionsSummary, Symbol("statusesSummary"), statusesSummary)
        setfield!(o, Symbol("statusesSummary"), statusesSummary)
        o
    end
end # type VirtualMachineScaleSetVMExtensionsSummary

const _property_map_VirtualMachineScaleSetVMExtensionsSummary = Dict{Symbol,Symbol}(Symbol("name")=>Symbol("name"), Symbol("statusesSummary")=>Symbol("statusesSummary"))
const _property_types_VirtualMachineScaleSetVMExtensionsSummary = Dict{Symbol,String}(Symbol("name")=>"String", Symbol("statusesSummary")=>"Vector{VirtualMachineStatusCodeCount}")
Base.propertynames(::Type{ VirtualMachineScaleSetVMExtensionsSummary }) = collect(keys(_property_map_VirtualMachineScaleSetVMExtensionsSummary))
Swagger.property_type(::Type{ VirtualMachineScaleSetVMExtensionsSummary }, name::Symbol) = Union{Nothing,eval(Meta.parse(_property_types_VirtualMachineScaleSetVMExtensionsSummary[name]))}
Swagger.field_name(::Type{ VirtualMachineScaleSetVMExtensionsSummary }, property_name::Symbol) =  _property_map_VirtualMachineScaleSetVMExtensionsSummary[property_name]

function check_required(o::VirtualMachineScaleSetVMExtensionsSummary)
    true
end

function validate_property(::Type{ VirtualMachineScaleSetVMExtensionsSummary }, name::Symbol, val)
end


text:
struct MultiChannelLayout
    nobs_per_channel :: Array{Int, 1}
    nchannels :: Int
    nobs :: Int
    # Inner constructor to validate fields
    MultiChannelLayout(nobs_per_channel) =
        any([p <= 0 for p in nobs_per_channel]) ? throw(ArgumentError("Must have non-empty channels")) :
             new(nobs_per_channel, size(nobs_per_channel)[1], sum(nobs_per_channel))
end

struct MultiChannelData
    values :: Array{Array{Array{Float64, 1}, 1}, 1}
    channellayout :: MultiChannelLayout
    nsamples :: Int
    obs_ids :: Array{Int, 1}
    sample_ids :: Array{Int, 1}
    function check_conds(values, channellayout, nsamples; obs_ids=[], sample_ids=[])
        if nsamples <= 0
            throw(ArgumentError("Must have non-zero samples"))
        elseif size(values)[1] != nsamples
            throw(ArgumentError("Values don't match number of samples"))
        elseif any([size(samp)[1] != channellayout.nchannels for samp in values])
            throw(ArgumentError("Values don't match number of channels"))
        elseif any([any([size(c)[1] != nchan for (c, nchan) in zip(samp, channellayout.nobs_per_channel)]) for samp in values])
            throw(ArgumentError("Values don't match number of observations per channel"))
        else
            if length(obs_ids) == 0 && length(sample_ids) == 0
                new(values, channellayout, nsamples, [], [])
            else
                new(values, channellayout, nsamples, obs_ids, sample_ids)
            end
        end
    end
    #MultiChannelData(values, channellayout, nsamples) = check_conds(values, channellayout, nsamples)
    MultiChannelData(values, channellayout, nsamples; obs_ids=[], sample_ids=[]) = check_conds(
        values, channellayout, nsamples; obs_ids, sample_ids)
end

struct MultiChannelFactorLayout
    nchannelshared :: Int # Number of channel-shared factors
    nspecific_per_channel :: Array{Int, 1}  # Number of channel-specific factors for each chan.
    nchannels :: Int # Number of channels
    nchannelspecific :: Int # Total channel-specific factors
    nfactors :: Int # Total number of factors
    MultiChannelFactorLayout(nchannelshared, nspecific_per_channel) = new(
        nchannelshared,
        nspecific_per_channel,
        size(nspecific_per_channel)[1],
        sum(nspecific_per_channel),
        nchannelshared + sum(nspecific_per_channel),
        )
end

struct MultiChannelFactors
    factorlayout :: MultiChannelFactorLayout
    fcommon :: Matrix{Float64}
    funique :: Matrix{Float64}
    funique_by_channel :: Array{Matrix{Float64}, 1}
    factors :: Matrix{Float64}
    sample_ids :: Array{Int, 1}
    function check_conds(factorlayout, fcommon, funique; sample_ids)
        funique_by_channel = extract_horizontal_blocks(funique, factorlayout.nspecific_per_channel; get_view=false)
        factors = vcat(fcommon, funique)
        if size(fcommon)[1] != factorlayout.nchannelshared
            throw(ArgumentError("Number of laid-out common factors does not match supplied"))
        elseif size(funique)[1] != factorlayout.nchannelspecific
            throw(ArgumentError("Number of laid-out unique factors does not match supplied"))
        elseif size(fcommon)[2] != size(funique)[2]
            throw(ArgumentError("Length of common and unique factors do not match"))
        else
            new(
            factorlayout,
            fcommon,
            funique,
            funique_by_channel,
            factors,
            sample_ids
            )
        end
    end
    #MultiChannelFactors(factorlayout, fcommon, funique) = check_conds(factorlayout, fcommon, funique)
    MultiChannelFactors(factorlayout, fcommon, funique; sample_ids=[]) = check_conds(factorlayout, fcommon, funique; sample_ids)
end

struct MCFMParams
    factorlayout :: MultiChannelFactorLayout
    channellayout :: MultiChannelLayout
    H :: Matrix{Float64}
    G :: Matrix{Float64}
    Σ :: Matrix{Float64}
    H_by_channel :: Array{Matrix{Float64},1}
    G_by_channel :: Array{Matrix{Float64},1}
    Σ_by_channel :: Array{Matrix{Float64},1}
    loading :: Matrix{Float64}
    obs_ids :: Array{Int, 1}
    function check_conds(factorlayout, channellayout, H, G, Σ; obs_ids)
        #TODO: Validate shapes of H/G/Sig by factorlayout
        H_by_channel = extract_horizontal_blocks(H, channellayout.nobs_per_channel; get_view=false)
        G_by_channel = extract_diagonal_blocks(G, channellayout.nobs_per_channel, factorlayout.nspecific_per_channel; get_view=false)
        Σ_by_channel = extract_diagonal_blocks(Σ, channellayout.nobs_per_channel, channellayout.nobs_per_channel; get_view=false)
        loading = hcat(H, G)
        new(factorlayout, channellayout, H, G, Σ, H_by_channel, G_by_channel, Σ_by_channel, loading, obs_ids)
    end
    MCFMParams(factorlayout, channellayout, H, G, Σ; obs_ids=[]) = check_conds(factorlayout, channellayout, H, G, Σ; obs_ids)
end

struct MCFMHistory
    niter :: Int
    maxiter :: Int
    tol :: Float64
    parampath :: Array{MCFMParams, 1}
end

struct MCFMFit
    finalparams :: MCFMParams
    history :: MCFMHistory
end


text:
struct GitLink
    # Git stuff
    root_dir::String
    remote_url::String

    # System data
    dat::Dict{Any, Any}

    GitLink(root_dir, remote_url) = new(string(root_dir), string(remote_url), Dict{Any, Any}())
    GitLink(;root_dir, remote_url) = GitLink(root_dir, remote_url)
end

# overwrite base
import Base.show

function Base.show(io::IO, gl::GitLink)
    println(io, "GitLink(;")
    println(io, "   root_dir = \"", gl.root_dir, "\",")
    println(io, "   remote_url = \"", gl.remote_url, "\"")
    println(io, ")")
end

import Base.get!
get!(f::Function, gl::GitLink, key) = get!(f, gl.dat, key)
get!(gl::GitLink, key, val) = get!(gl.dat, key, val)

import Base.get
get(f::Function, gl::GitLink, key) = get(f, gl.dat, key)
get(gl::GitLink, key, val) = get(gl.dat, key, val)

set!(gl::GitLink, key, val) = (gl.dat[key] = val)
set!(f::Function, gl::GitLink, key) = (gl.dat[key] = f())

remote_url(gl::GitLink) = gl.remote_url


text:
mutable struct Model
    body::Expr
    params::Expr
    arguments::Tuple
    conditionings::Dict{Symbol,Any}
    observed::Set{Symbol}

    Model(body::Expr) = new(body, :(), (), Dict(), Set())
    Model(body::Expr, params::Expr) = new(body, params, (), Dict(), Set())
    Model(body::Expr, params::Expr, arguments::Tuple, conditionings::Dict{Symbol,Any}, observed::Set{Symbol}) = new(body, params, arguments, conditionings, observed)
end

macro model(body::Expr)
    Model(prettify(body))
end

macro model(param::Symbol, body::Expr)
    Model(prettify(body), :(($param,)))
end

macro model(params::Expr, body::Expr)
    Model(prettify(body), params)
end

function copy(model::Model)
    Model(model.body, model.params, model.arguments, Dict(model.conditionings), Set(model.observed))
end

function condition{T}(conds::Vararg{Pair{Symbol,T}})
    (model::Model) -> let
        newmodel = copy(model)
        for cond in conds
            newmodel.conditionings[cond.first] = cond.second
        end
        newmodel
    end
end

function observe(obs::Vararg{Symbol})
    (model::Model) -> let
        newmodel = copy(model)
        union!(newmodel.observed, Set(obs))
        newmodel
    end
end

function (model::Model)(args :: Vararg)
    newmodel = copy(model)
    newmodel.arguments = args
    newmodel
end


text:
mutable struct Asset
    copyright::Union{Nothing,String}
    generator::Union{Nothing,String}
    version::String
    minVersion::Union{Nothing,String}
    extensions::Union{Nothing,Dict}
    extras::Union{Nothing,Dict}
    function Asset()
        obj = new()
        obj.copyright = nothing
        obj.generator = nothing
        obj.minVersion = nothing
        obj.extensions = nothing
        obj.extras = nothing
        obj
    end
end

JSON3.StructType(::Type{Asset}) = JSON3.Mutable()
StructTypes.omitempties(::Type{Asset}) = (:copyright, :generator, :minVersion, :extensions, :extras)


text:
module AeroAcoustics
using Distances
using LinearAlgebra
using Statistics
using NLsolve
using Parameters
using LazyArrays
import DSP

import Base.length,
       Base.push!,
       Base.reshape

export Environment,
       FreqArray,
       SPL,
       steeringvectors,
       steeringvectors!,
       csm,
       beamforming,
       psf,
       damas,
       cleanSC,
       sourceintegration,
       octavebandlimits,
       octavebands,
       narrow2oct,
       enbw,
       coherence_weights,
       DR!,
       DR

include("utils.jl")
include("types.jl")
include("steeringvectors.jl")
include("csm.jl")
include("beamforming.jl")
include("psf.jl")
include("damas.jl")
include("cleansc.jl")
include("sourceintegration.jl")

reshape(x::FreqArray,dims::Vararg{Int64,N}) where N = FreqArray(reshape(x.arr,dims...),x.fc)

"""
    narrow2oct(x::FreqArray,n,nomial::Bool=true;psd=false)

Sum narrow band spectrum to 1/n octave bands given narrow band frequencies `f` in x.
"""
function narrow2oct(x::FreqArray,n;nomial::Bool=true,psd::Bool=false)
    fc = octavebands(n,extrema(x.fc);nomial=nomial)
    fl, fu = octavebandlimits(fc,n)
    out = Array{Union{Missing, Float64},2}(undef,size(x,1),length(fc))
    for i in 1:length(fc)
        inds = selectfreqs(x.fc,(fl[i],fu[i]))
        if any(inds .== true)
            if psd
                out[:,i] = sum(x[:,inds],dims=2)/(x.fc[2]-x.fc[1])
            else
                out[:,i] = sum(x[:,inds],dims=2)
            end
        else
            @warn "No data in frequency band $(fc[i]) replacing with missing"
            out[:,i] = Array{Missing,1}(undef,size(x,1))
        end
    end
    return FreqArray(out,fc)
end

"""
    DR!(CSM::FreqArray)

Apply diagonal removal to CSM.
"""
function DR!(CSM::FreqArray)
    T = eltype(CSM)
    M,M,Nf = size(CSM)
    for j = 1:Nf
        for i = 1:M
            CSM.arr[i,i,j] = zero(T)
        end
    end
end

"""
    DR(CSM::FreqArray)

Apply diagonal removal to CSM and new array.
"""
function DR(CSM::FreqArray)
    T = eltype(CSM)
    M,M,Nf = size(CSM)
    CSMd = copy(CSM)
    for j = 1:Nf
        for i = 1:M
            CSMd[i,i,j] = zero(T)
        end
    end
    return FreqArray(CSMd,CSM.fc)
end

end


text:
type Sigmoid{A<:AbstractArray} <: Pointwise{A}
    output :: A
    gradInput :: A
end

Sigmoid(T::Type, s...) = Sigmoid(Array(T, s...), Array(T, s...))

sigmoid(x) = 1 / (1 + exp(-x))
sigmoidPrime(x) = sigmoid(x) * (1 - sigmoid(x))

function updateOutput!{A}(self :: Sigmoid{A}, input :: A)
    map!(sigmoid, self.output, input)
    return self.output
end

function updateGradInput!{A}(self :: Sigmoid{A}, input :: A, gradOutput :: A)
    map!(sigmoidPrime, self.gradInput, gradOutput)
    return self.gradInput
end

export Sigmoid


text:
function is_stationary(v, th, w)
    length(v) < w && return false
    @views pv = v[end - w + 1:end]
    m = Statistics.mean(pv)
    m = iszero(m) ? one(m) : m
    s = Statistics.std(pv)
    return s <= abs(m) * th
end


text:
import Base.zero

export
    AV_PERM_READ,
    AV_PERM_WRITE,
    AV_PERM_PRESERVE,
    AV_PERM_REUSE,
    AV_PERM_REUSE2,
    AV_PERM_NEG_LINESIZES,
    AV_PERM_ALIGN,
    AVFILTER_ALIGN,
    AVFILTER_FLAG_DYNAMIC_INPUTS,
    AVFILTER_FLAG_DYNAMIC_OUTPUTS,
    AVFILTER_FLAG_SLICE_THREADS,
    AVFILTER_FLAG_SUPPORT_TIMELINE_GENERIC,
    AVFILTER_FLAG_SUPPORT_TIMELINE_INTERNAL,
    AVFILTER_FLAG_SUPPORT_TIMELINE,
    AVFILTER_THREAD_SLICE,
    AVFILTER_CMD_FLAG_ONE,
    AVFILTER_CMD_FLAG_FAST,
    AVFilterContext,
    AVFilterLink,
    AVFilterPad,
    AVFilterFormats,
    Array_8_Ptr,
    Array_8_Cint,
    AVFilterBuffer,
    AVFilterBufferRefAudioProps,
    AVFilterBufferRefVideoProps,
    AVFilterBufferRef,
    AVFilter,
    AVFilterInternal,
    AVFilterGraphInternal,
    AVFilterGraph,
    ANONYMOUS_1,
    AVFILTER_AUTO_CONVERT_ALL,
    AVFILTER_AUTO_CONVERT_NONE,
    AVFilterInOut,
    AV_BUFFERSINK_FLAG_PEEK,
    AV_BUFFERSINK_FLAG_NO_REQUEST,
    AVBufferSinkParams,
    AVABufferSinkParams,
    ANONYMOUS_2,
    AV_BUFFERSRC_FLAG_NO_CHECK_FORMAT,
    AV_BUFFERSRC_FLAG_NO_COPY,
    AV_BUFFERSRC_FLAG_PUSH,
    AV_BUFFERSRC_FLAG_KEEP_REF

const AV_PERM_READ = 0x01
const AV_PERM_WRITE = 0x02
const AV_PERM_PRESERVE = 0x04
const AV_PERM_REUSE = 0x08
const AV_PERM_REUSE2 = 0x10
const AV_PERM_NEG_LINESIZES = 0x20
const AV_PERM_ALIGN = 0x40
const AVFILTER_ALIGN = 16
const AVFILTER_FLAG_DYNAMIC_INPUTS = 1 << 0
const AVFILTER_FLAG_DYNAMIC_OUTPUTS = 1 << 1
const AVFILTER_FLAG_SLICE_THREADS = 1 << 2
const AVFILTER_FLAG_SUPPORT_TIMELINE_GENERIC = 1 << 16
const AVFILTER_FLAG_SUPPORT_TIMELINE_INTERNAL = 1 << 17
const AVFILTER_FLAG_SUPPORT_TIMELINE = AVFILTER_FLAG_SUPPORT_TIMELINE_GENERIC | AVFILTER_FLAG_SUPPORT_TIMELINE_INTERNAL
const AVFILTER_THREAD_SLICE = 1 << 0
const AVFILTER_CMD_FLAG_ONE = 1
const AVFILTER_CMD_FLAG_FAST = 2

typealias AVFilterContext AVFilterContext
typealias AVFilterLink AVFilterLink
typealias AVFilterPad AVFilterPad
typealias AVFilterFormats Void

immutable Array_8_Ptr
    d1::Ptr{UInt8}
    d2::Ptr{UInt8}
    d3::Ptr{UInt8}
    d4::Ptr{UInt8}
    d5::Ptr{UInt8}
    d6::Ptr{UInt8}
    d7::Ptr{UInt8}
    d8::Ptr{UInt8}
end

zero(::Type{Array_8_Ptr}) = Array_8_Ptr(fill(C_NULL,8)...)

immutable Array_8_Cint
    d1::Cint
    d2::Cint
    d3::Cint
    d4::Cint
    d5::Cint
    d6::Cint
    d7::Cint
    d8::Cint
end

zero(::Type{Array_8_Cint}) = Array_8_Cint(fill(zero(Cint),8)...)

immutable AVFilterBuffer
    data::Array_8_Ptr
    extended_data::Ptr{Ptr{UInt8}}
    linesize::Array_8_Cint
    priv::Ptr{Void}
    free::Ptr{Void}
    format::Cint
    w::Cint
    h::Cint
    refcount::UInt32
end

immutable AVFilterBufferRefAudioProps
    channel_layout::UInt64
    nb_samples::Cint
    sample_rate::Cint
    channels::Cint
end

immutable AVFilterBufferRefVideoProps
    w::Cint
    h::Cint
    sample_aspect_ratio::AVRational
    interlaced::Cint
    top_field_first::Cint
    pict_type::AVPictureType
    key_frame::Cint
    qp_table_linesize::Cint
    qp_table_size::Cint
    qp_table::Ptr{Int8}
end

immutable AVFilterBufferRef
    buf::Ptr{AVFilterBuffer}
    data::Array_8_Ptr
    extended_data::Ptr{Ptr{UInt8}}
    linesize::Array_8_Cint
    video::Ptr{AVFilterBufferRefVideoProps}
    audio::Ptr{AVFilterBufferRefAudioProps}
    pts::Int64
    pos::Int64
    format::Cint
    perms::Cint
    _type::AVMediaType
    metadata::Ptr{AVDictionary}
end

immutable AVFilter
    name::Ptr{UInt8}
    description::Ptr{UInt8}
    inputs::Ptr{AVFilterPad}
    outputs::Ptr{AVFilterPad}
    priv_class::Ptr{AVClass}
    flags::Cint
    init::Ptr{Void}
    init_dict::Ptr{Void}
    uninit::Ptr{Void}
    query_formats::Ptr{Void}
    priv_size::Cint
    next::Ptr{AVFilter}
    process_command::Ptr{Void}
    init_opaque::Ptr{Void}
end

typealias AVFilterInternal Void
typealias AVFilterGraphInternal Void

immutable AVFilterGraph
    av_class::Ptr{AVClass}
    filter_count_unused::UInt32
    filters::Ptr{Ptr{AVFilterContext}}
    scale_sws_opts::Ptr{UInt8}
    resample_lavr_opts::Ptr{UInt8}
    nb_filters::UInt32
    thread_type::Cint
    nb_threads::Cint
    internal::Ptr{AVFilterGraphInternal}
    opaque::Ptr{Void}
    execute::Ptr{avfilter_execute_func}
    aresample_swr_opts::Ptr{UInt8}
    sink_links::Ptr{Ptr{AVFilterLink}}
    sink_links_count::Cint
    disable_auto_convert::UInt32
end

# begin enum ANONYMOUS_1
typealias ANONYMOUS_1 Cint
const AVFILTER_AUTO_CONVERT_ALL = Int32(0)
const AVFILTER_AUTO_CONVERT_NONE = Int32(-1)
# end enum ANONYMOUS_1

immutable AVFilterInOut
    name::Ptr{UInt8}
    filter_ctx::Ptr{AVFilterContext}
    pad_idx::Cint
    next::Ptr{AVFilterInOut}
end

const AV_BUFFERSINK_FLAG_PEEK = 1
const AV_BUFFERSINK_FLAG_NO_REQUEST = 2

immutable AVBufferSinkParams
    pixel_fmts::Ptr{AVPixelFormat}
end

immutable AVABufferSinkParams
    sample_fmts::Ptr{AVSampleFormat}
    channel_layouts::Ptr{Int64}
    channel_counts::Ptr{Cint}
    all_channel_counts::Cint
    sample_rates::Ptr{Cint}
end

# begin enum ANONYMOUS_2
typealias ANONYMOUS_2 UInt32
const AV_BUFFERSRC_FLAG_NO_CHECK_FORMAT = UInt32(1)
const AV_BUFFERSRC_FLAG_NO_COPY = UInt32(2)
const AV_BUFFERSRC_FLAG_PUSH = UInt32(4)
const AV_BUFFERSRC_FLAG_KEEP_REF = UInt32(8)
# end enum ANONYMOUS_2


text:
using ITensors
using Test

# Wrap an ITensor with pairs of primed and
# unprimed indices to pass to davidson
struct ITensorMap
  A::ITensor
end
Base.eltype(M::ITensorMap) = eltype(M.A)
Base.size(M::ITensorMap) = dim(findinds(M.A,("",0)))
(M::ITensorMap)(v::ITensor) = noprime(M.A*v)

@testset "Complex davidson" begin
  d = 10
  i = Index(d,"i")
  A = randomITensor(Complex,i,prime(i))
  A = mapprime(A*mapprime(dag(A),0,2),2,1)
  M = ITensorMap(A)

  v = randomITensor(i)
  λ,v = davidson(M,v;maxiter=10)
  @test M(v)≈λ*v

  v = randomITensor(Complex, i)
  λ,v = davidson(M,v;maxiter=10)
  @test M(v)≈λ*v

end


text:
using Test, JLLPrefixes, Base.BinaryPlatforms, Pkg
using JLLPrefixes: PkgSpec

@testset "JLL collection" begin
    function check_zstd_jll(zstd_pkgspec, zstd_artifacts)
        # Ensure this pkgspec is named Zstd_jll
        @test zstd_pkgspec.name == "Zstd_jll"

        # It only had one artifact to download, and it exists
        @test length(zstd_artifacts) == 1
        @test isdir(zstd_artifacts[1])
        @test isfile(joinpath(zstd_artifacts[1], "include", "zstd.h"))
    end

    @testset "Zstd_jll (native)" begin
        # Start with a simple JLL with no dependencies
        artifact_paths = collect_artifact_paths(["Zstd_jll"]; verbose=true)

        # There was only one JLL downloaded, and it was Zstd_jll
        @test length(artifact_paths) == 1
        zstd_pkgspec, zstd_artifacts = first(artifact_paths)
        check_zstd_jll(zstd_pkgspec, zstd_artifacts)
    end

    # Do another simple JLL installation, but this time for a few different architectures
    for platform in [Platform("aarch64", "linux"), Platform("x86_64", "macos"), Platform("i686", "windows")]
        @testset "Zstd_jll ($(platform))" begin
            artifact_paths = collect_artifact_paths(["Zstd_jll"]; platform)
            check_zstd_jll(first(artifact_paths)...)

            # Test that we're getting the kind of dynamic library we expect
            artifact_dir = first(first(values(artifact_paths)))
            if os(platform) == "windows"
                libdir = "bin"
                libname = "libzstd-1.dll"
            elseif os(platform) == "macos"
                libdir = "lib"
                libname = "libzstd.1.dylib"
            else
                libdir = "lib"
                libname = "libzstd.so.1"
            end
            @test isfile(joinpath(artifact_dir, libdir, libname))
        end
    end

    # Test that we can request a particular version of Zstd_jll
    linux64 = Platform("x86_64", "linux")
    @testset "Zstd_jll ($(linux64), v1.4.2)" begin
        artifact_paths = collect_artifact_paths([PkgSpec(;name="Zstd_jll", version=v"1.4.2")]; platform=linux64)

        # There was only one JLL downloaded, and it was Zstd_jll
        @test length(artifact_paths) == 1
        zstd_pkgspec, zstd_artifacts = first(artifact_paths)
        check_zstd_jll(zstd_pkgspec, zstd_artifacts)

        # Ensure that this is actually version 1.4.2
        artifact_dir = first(first(values(artifact_paths)))
        @test isfile(joinpath(artifact_dir, "lib", "libzstd.so.1.4.2"))
    end

    # Kick it up a notch; start involving dependencies
    @testset "XML2_jll ($(linux64), dependencies)" begin
        # Lock XML2_jll to v2.9 in case it adds more dependencies in the future
        artifact_paths = collect_artifact_paths([PkgSpec(;name="XML2_jll", version=v"2.9.12")]; platform=linux64)

        @test length(artifact_paths) == 3
        @test sort([p.name for p in keys(artifact_paths)]) == ["Libiconv_jll", "XML2_jll", "Zlib_jll"]
    end

    # Install two packages that have nothing to do with eachother at the same time
    @testset "Bzip2_jll + Zstd_jll" begin
        artifact_paths = collect_artifact_paths(["Bzip2_jll", "Zstd_jll"])
        @test length(artifact_paths) == 2
        @test sort([p.name for p in keys(artifact_paths)]) == ["Bzip2_jll", "Zstd_jll"]
    end

    # Test stdlibs across versions
    for (GMP_soversion, julia_version) in (("10.3.2", v"1.5"), ("10.4.0", v"1.6"), ("10.4.1", v"1.7"))
        @testset "GMP_jll (Julia $(julia_version))" begin
            artifact_paths = collect_artifact_paths(["GMP_jll"]; platform=Platform("x86_64", "linux"; julia_version))
            @test length(artifact_paths) == 1
            gmp_artifact_dir = only(first(values(artifact_paths)))
            @test isfile(joinpath(gmp_artifact_dir, "lib", "libgmp.so.$(GMP_soversion)"))
        end
    end

    # Test "impossible" situations via `julia_version == nothing`
    @testset "Impossible Constraints" begin
        # We can't naively install GMP v6.1.2 and MPFR v4.1.0, because those are
        # from conflicting Julia versions, and the Pkg resolver doesn't like that
        for julia_version in (v"1.5", v"1.6")
            @test_throws Pkg.Resolve.ResolverError collect_artifact_paths([
                PkgSpec(;name="GMP_jll",  version=v"6.1.2"),
                PkgSpec(;name="MPFR_jll", version=v"4.1.0"),
            ]; platform=Platform("x86_64", "linux"; julia_version))
        end

        # So we must pass julia_version == nothing, as is the case in our `linux64` object
        artifact_paths = collect_artifact_paths([
            PkgSpec(;name="GMP_jll",  version=v"6.1.2"),
            PkgSpec(;name="MPFR_jll", version=v"4.1.0"),
        ]; platform=linux64)
        @test length(artifact_paths) == 2
        @test sort([p.name for p in keys(artifact_paths)]) == ["GMP_jll", "MPFR_jll"]
    end

    # Test adding something that doesn't exist on a certain platform
    @testset "Platform Incompatibility" begin
        @test_logs (:warn, r"Dependency Libuuid_jll does not have a mapping for artifact Libuuid for platform x86_64-apple-darwin") begin
            artifact_paths = collect_artifact_paths(["Libuuid_jll"]; platform=Platform("x86_64", "macos"), verbose=true)
            @test isempty(artifact_paths)
        end
    end
end

exe = ""
if Sys.iswindows()
    exe = ".exe"
end

@testset "Symlinking" begin
    # Grab a big JLL, like FFMPEG, symlink them all together, and ensure that all the files exist
    @testset "FFMPEG" begin
        artifact_paths = collect_artifact_paths(["FFMPEG_jll"])
        mktempdir() do prefix
            symlink_artifact_paths(prefix, artifact_paths)

            # Ensure that a bunch of tools we expect to be installed are, in fact, installed
            for tool in ("ffmpeg", "bzcat", "fc-cache", "iconv", "x264", "x265", "xslt-config")
                @test isfile(joinpath(prefix, "bin", "$(tool)$(exe)"))
                @test isfile(abspath(joinpath(prefix, "bin", "$(tool)$(exe)")))
            end
        end
    end
end

@testset "Copying" begin
    # Grab a big JLL, like FFMPEG, symlink them all together, and ensure that all the files exist
    # Also, this should allow us to actually run `ffmpeg`!
    @testset "FFMPEG" begin
        artifact_paths = collect_artifact_paths(["FFMPEG_jll"])
        mktempdir() do prefix
            copy_artifact_paths(prefix, artifact_paths)

            # Ensure that a bunch of tools we expect to be installed are, in fact, installed
            for tool in ("ffmpeg", "bzcat", "fc-cache", "iconv", "x264", "x265", "xslt-config")
                # Use `eval` so that the test failure shows which tools fail
                tool_name = string(tool, exe)
                @eval @test isfile(joinpath($(prefix), "bin", $(tool_name)))
            end

            run(`$(joinpath(prefix, "bin", "ffmpeg$(exe)")) -version`)
        end
    end
end


text:
module Pollen

using HTTP, Gumbo, CSV, Cascadia, DataFrames, JSONTables, Dates, Query
using Pipe:@pipe
using Cascadia:matchFirst

const CSV_POLLEN = joinpath(@__DIR__, "levels.csv")
const JSON_POLLEN = joinpath(@__DIR__, "levels.json")

function fetch()::NamedTuple
    url = ENV["POLLEN_API_URL"]
    response = HTTP.get(url)
    html = parsehtml(String(response))
    dates = parsedate(html.root)
    elm, willow, poplar, hazel, alder, oak, beech, birch, mugwort, ragweed, plantain, sorrel, rye, grass = parsetable(html.root)

    (elm = elm,
    willow = willow,
    poplar = poplar,
    hazel = hazel,
    alder = alder,
    oak = oak,
    beech = beech,
    birch = birch,
    mugwort = mugwort,
    ragweed = ragweed,
    plantain = plantain,
    sorrel = sorrel,
    rye = rye,
    grass = grass,
    date = dates)
end

function record()
    persisted = CSV.read(CSV_POLLEN, DataFrame)
    df = @pipe fetch() |> DataFrame |> vcat(_, persisted) |> unique(_, :date) |> sort!(_, :date, rev=true)
    df |> CSV.write(CSV_POLLEN)

    open(JSON_POLLEN, "w") do f
        write(f, arraytable(df |> @mutate(formatted_date = format_date(_.date)) |> DataFrame))
    end
end

format_date(date::Date)::String = "$(Dates.dayname(date)), $(Dates.monthname(date)) $(Dates.day(date))"

function parsedate(root)::Vector{Date}
    datestring = matchFirst(sel"#date0 p", root)[1].text
    parts = match(r"(\d+)\.(\d+)\.", datestring).captures
    current = Date(Dates.year(today()), parse(Int, parts[2]), parse(Int, parts[1]))
    collect(current:Dates.Day(1):current + Dates.Day(6))
end

function parsetable(root)
    table = matchFirst(sel"#pollen_tabelle", root)
    parsecategories(table, ("elm_text",
                            "wil_text",
                            "pop_text",
                            "haz_text",
                            "ald_text",
                            "oak_text",
                            "bee_text",
                            "bir_text",
                            "mug_text",
                            "hog_text",
                            "pla_text",
                            "rum_text",
                            "rye_text",
                            "gra_text",))
end

parsecategories(table, selectors) = map(s -> parsecategory(table, s), selectors)

function parsecategory(table, selector)
    category = matchFirst(Selector("#$selector + .grid_item_pollen_tabelle"), table)
    map(b -> grade(b), dayvalues(category))
end

function dayvalues(category)
    [category.attributes["data-day0"],
     category.attributes["data-day1"],
     category.attributes["data-day2"],
     category.attributes["data-day3"],
     category.attributes["data-day4"],
     category.attributes["data-day5"],
     category.attributes["data-day6"]]
end

function grade(burden)
    burdens = ("noburden", "weakburden", "moderateburden", "strongburden")
    findfirst(b -> b == burden, burdens) - 1
end

end # module


text:
# # # RAMI Benchmarking
#  The goal of the RAMI4PILPS experiment[^1] is to evaluate different approaches by which Land Surface Models in larger Earth System Models quantify the radiation transfer within vegetation canopies. The RAMI4PILPS can be interpreted as a quality control mechanism used to:
#
# 1) quantify the errors in the radiative transfer scheme;
#
# 2) identify the impact that structural and spectral sub-grid variability may have on radiative transfer; and
#
# 3) verify the conservation of energy at the level of the surface, as well as inconsistencies arising from different levels of assumptions/simplifications.
#
# This approach involves direct comparison with reference solutions obtained from highly accurate 3D models identified during the third phase of the RAMI benchmarking exercise[^2].
#
# A set of 3D experiments compares the partitioning of incident solar energy into an absorbed (A) flux, a transmitted (T) flux component and the surface reflectance (R). The overall canopy structure for these test cases is reminiscent of open forest canopies with randomly oriented foliage, confined to spherical volumes located at varying heights above the ground.
# In here, we use a total of 36 test cases including various canopy density, soil brightness, and illumination conditions for the visible (VIS) (400-700nm) and near infra-red (NIR) (700-3000nm) spectral ranges.
#
#  [^1]: J.L. Widlowski, B. Pinty, M. Clerici, Y. Dai, M. De Kauwe, K. de Ridder, A. Kallel, H. Kobayashi, T. Lavergne, W. Ni-Meister, A. Olchev, T. Quaife, S. Wang, W. Yang, Y. Yang, and H. Yuan (2011), RAMI4PILPS: An intercomparison of formulations for the partitioning of solar radiation in land surface models, Journal of Geophysical Research, 116, G02019, 25, DOI: 10.1029/2010JG001511.
#
#  [^2]: Widlowski, J-L., M. Taberner, B. Pinty, V. Bruniquel-Pinel, M. Disney, R. Fernandes, J.-P. Gastellu-Etchegorry, N. Gobron, A. Kuusk, T. Lavergne, S. Leblanc, P. Lewis, E. Martin, M. Mottus, P. J. R. North, W. Qin, M.Robustelli, N. Rochdi, R.Ruiloba, C.Soler, R.Thompson, W. Verhoef, M. M.Verstraete, and D. Xie (2007), 'The third RAdiation transfer Model Intercomparison (RAMI) exercise: Documenting progress in canopy reflectance modelling', Journal of Geophysical Research, 112, D09111, 28, DOI: 10.1029/2006JD007821.
#

#### Use Julia Plots package and switch to plotly js option:
using PyPlot

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

## First, we include Revise (good for debugging) and Parameters (tools for structures)

##using Revise
using Parameters
##----------------------------------------------------------------------------

## Now include the Land modules

using Land
using Land.CanopyRT
##----------------------------------------------------------------------------

##Defining all reference values for the Sparse case

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

const FT = Float32

wl_set = create_wl_para_set(FT)
leaf = create_leaf_bio(FT, wl_set.nwl, wl_set.nWlE, wl_set.nWlF);
canopy_rt = Canopy4RT{FT, 20, 3.0}()
canRad_rt = CanopyRadiation{FT, wl_set.nwl, wl_set.nWlF, length(canopy_rt.litab), length(canopy_rt.lazitab), canopy_rt.nlayers}()
canOpt_rt = create_canopy_optical(FT, wl_set.nwl, canopy_rt.nlayers, length(canopy_rt.lazitab), length(canopy_rt.litab); using_marray=false)
sunRad_rt = create_incoming_radiation(FT, wl_set.swl);

soil = SoilOpti{FT}(wl_set.wl, FT(0.2)*ones(FT, length(wl_set.wl)), FT[0.1], FT(290.0))
angles = SolarAngles{FT}()

arrayOfLeaves = [create_leaf_bio(FT, wl_set.nwl, wl_set.nWlE, wl_set.nWlF) for i in 1:canopy_rt.nlayers]
for i in 1:canopy_rt.nlayers
    fluspect!(arrayOfLeaves[i], wl_set)
end

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

RAMI_SZA = [27.,60.,83.]

RAMI_fabsRed_050_BLK =  [0.09380509999999999, 0.16259713, 0.53931207]
RAMI_frefRed_050_BLK =  [0.00330673, 0.00517598, 0.01626682]
RAMI_ftranRed_050_BLK =  [0.90288817, 0.83222689, 0.44442110999999995]

RAMI_fabsRed_050_MED =  [0.10897124, 0.17760124000000002, 0.54764719]
RAMI_frefRed_050_MED =  [0.09759354, 0.09107608, 0.06177913]
RAMI_ftranRed_050_MED =  [0.90337609, 0.83265704, 0.44469279]

RAMI_fabsRed_050_SNW =  [0.21471034, 0.28200132, 0.60564705]
RAMI_frefRed_050_SNW =  [0.7526521700000001, 0.6879087300000001, 0.37825442000000004]
RAMI_ftranRed_050_SNW =  [0.90659694, 0.83583194, 0.44718138]
#----------------------------------------------------------------------------

function RAMI_case(LAI, soil_albedo, clumping_index)

  soil.albedo_SW[:] .=soil_albedo;
  ##Clumping index
  canopy_rt.Ω = 1.0
  ##Viewing Zenith Angle in degrees
  angles.tto=0.0
  ##Leaf Area index
  canopy_rt.LAI=LAI

  reflRed_SZA = []
  absRed_SZA = []
  transRed_SZA = []

  for SZA=0.0:1:85
    angles.tts=SZA

    fluspect!(leaf, wl_set);
    compute_canopy_geometry!(canopy_rt, angles, canOpt_rt)
    compute_canopy_matrices!(arrayOfLeaves, canOpt_rt);

   ## leaf reflectance RED
   leaf.ρ_SW[28] = 0.0735
   ## leaf transmittance
   leaf.τ_SW[28]= 0.0566

    ##Setting all diffuse to zero
    sunRad_rt.E_diffuse[28] = 0.0

    simulate_short_wave!(canopy_rt, canOpt_rt, canRad_rt, sunRad_rt, soil);
    push!(reflRed_SZA, canRad_rt.alb_direct[28])
    push!(absRed_SZA, (sum(canRad_rt.netSW_shade,dims=2)[28,1].+sum(canRad_rt.netSW_sunlit,dims=2)[28,1])./(sunRad_rt.E_diffuse[28].+sunRad_rt.E_direct[28]))
    push!(transRed_SZA,  (canOpt_rt.Es_[28,end] .+ canRad_rt.E_down[28,end])./(sunRad_rt.E_diffuse[28].+sunRad_rt.E_direct[28]))
  end

  ######## Clumped Case

  reflRed_clump_SZA = []
  absRed_clump_SZA = []
  transRed_clump_SZA = []

  ##Clumping index
  canopy_rt.Ω = clumping_index

  for SZA=0.0:1:85
    angles.tts=SZA

    fluspect!(leaf, wl_set);
    compute_canopy_geometry!(canopy_rt, angles, canOpt_rt)
    compute_canopy_matrices!(arrayOfLeaves, canOpt_rt);

    simulate_short_wave!(canopy_rt, canOpt_rt, canRad_rt, sunRad_rt, soil);
    push!(reflRed_clump_SZA, canRad_rt.alb_direct[28])
    push!(absRed_clump_SZA, (sum(canRad_rt.netSW_shade,dims=2)[28,1].+sum(canRad_rt.netSW_sunlit,dims=2)[28,1])./(sunRad_rt.E_diffuse[28].+sunRad_rt.E_direct[28]))
    push!(transRed_clump_SZA,  (canOpt_rt.Es_[28,end] .+ canRad_rt.E_down[28,end])./(sunRad_rt.E_diffuse[28].+sunRad_rt.E_direct[28]))

  end

  return reflRed_SZA,absRed_SZA,transRed_SZA,reflRed_clump_SZA,absRed_clump_SZA,transRed_clump_SZA

end;

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

#TODO nest those into a loop!

##Sparse case with black soil
reflRed_SZA,absRed_SZA,transRed_SZA,reflRed_clump_SZA,absRed_clump_SZA,transRed_clump_SZA=RAMI_case(0.50265, 0.0, 0.365864235);

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

SZA=0:1:85

figure(figsize=(10,5))

subplot(1,2,1)
plot(SZA, reflRed_SZA,label=["reflectance"])
plot(SZA, absRed_SZA ,label=["absorptance"])
plot(SZA, transRed_SZA,label=["transmittance"])
scatter(RAMI_SZA, RAMI_frefRed_050_BLK)
scatter(RAMI_SZA, RAMI_fabsRed_050_BLK)
scatter(RAMI_SZA, RAMI_ftranRed_050_BLK)
title("050_BLK - Default")
ylabel("Radiation Partitioning")
xlabel("Sun Zenith Angle")
xlim(0.0, 90.)
ylim(-0.05, 1.0)
xticks(0:20:91.)
yticks(0:0.2:1.0)

subplot(1,2,2)
plot(SZA, reflRed_clump_SZA,label=["reflectance"])
plot(SZA, absRed_clump_SZA ,label=["absorptance"])
plot(SZA,transRed_clump_SZA,label=["transmittance"])
scatter(RAMI_SZA, RAMI_frefRed_050_BLK,label="RAMI reflectance")
scatter(RAMI_SZA, RAMI_fabsRed_050_BLK,label="RAMI absorptance")
scatter(RAMI_SZA, RAMI_ftranRed_050_BLK,label="RAMI transmittance")
title("Clumping")
ylabel("Radiation Partitioning")
xlabel("Sun Zenith Angle")
xlim(0.0, 90.)
ylim(-0.05, 1.0)
xticks(0:20:91.)
yticks(0:0.2:1.0)

legend()
gcf()

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

##Sparse case with medium soil
reflRed_SZA,absRed_SZA,transRed_SZA,reflRed_clump_SZA,absRed_clump_SZA,transRed_clump_SZA=RAMI_case(0.50265, 0.1217, 0.365864235);

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

SZA=0:1:85

figure(figsize=(10,5))

subplot(1,2,1)
plot(SZA, reflRed_SZA,label=["reflectance"])
plot(SZA, absRed_SZA ,label=["absorptance"])
plot(SZA, transRed_SZA,label=["transmittance"])
scatter(RAMI_SZA, RAMI_frefRed_050_MED)
scatter(RAMI_SZA, RAMI_fabsRed_050_MED)
scatter(RAMI_SZA, RAMI_ftranRed_050_MED)
title("050_MED - Default")
ylabel("Radiation Partitioning")
xlabel("Sun Zenith Angle")
xlim(0.0, 90.)
ylim(-0.05, 1.0)
xticks(0:20:91.)
yticks(0:0.2:1.0)

subplot(1,2,2)
plot(SZA, reflRed_clump_SZA,label=["reflectance"])
plot(SZA, absRed_clump_SZA ,label=["absorptance"])
plot(SZA,transRed_clump_SZA,label=["transmittance"])
scatter(RAMI_SZA, RAMI_frefRed_050_MED,label="RAMI reflectance")
scatter(RAMI_SZA, RAMI_fabsRed_050_MED,label="RAMI absorptance")
scatter(RAMI_SZA, RAMI_ftranRed_050_MED,label="RAMI transmittance")
title("Clumping")
ylabel("Radiation Partitioning")
xlabel("Sun Zenith Angle")
xlim(0.0, 90.)
ylim(-0.05, 1.0)
xticks(0:20:91.)
yticks(0:0.2:1.0)

legend()
gcf()

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

##Sparse case with snowy soil
reflRed_SZA,absRed_SZA,transRed_SZA,reflRed_clump_SZA,absRed_clump_SZA,transRed_clump_SZA=RAMI_case(0.50265, 0.9640, 0.365864235);

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

SZA=0:1:85

figure(figsize=(10,5))

subplot(1,2,1)
plot(SZA, reflRed_SZA,label=["reflectance"])
plot(SZA, absRed_SZA ,label=["absorptance"])
plot(SZA, transRed_SZA,label=["transmittance"])
scatter(RAMI_SZA, RAMI_frefRed_050_SNW)
scatter(RAMI_SZA, RAMI_fabsRed_050_SNW)
scatter(RAMI_SZA, RAMI_ftranRed_050_SNW)
title("050_SNW - Default")
ylabel("Radiation Partitioning")
xlabel("Sun Zenith Angle")
xlim(0.0, 90.)
ylim(-0.05, 1.0)
xticks(0:20:91.)
yticks(0:0.2:1.0)

subplot(1,2,2)
plot(SZA, reflRed_clump_SZA,label=["reflectance"])
plot(SZA, absRed_clump_SZA ,label=["absorptance"])
plot(SZA,transRed_clump_SZA,label=["transmittance"])
scatter(RAMI_SZA, RAMI_frefRed_050_SNW,label="RAMI reflectance")
scatter(RAMI_SZA, RAMI_fabsRed_050_SNW,label="RAMI absorptance")
scatter(RAMI_SZA, RAMI_ftranRed_050_SNW,label="RAMI transmittance")
title("Clumping")
ylabel("Radiation Partitioning")
xlabel("Sun Zenith Angle")
xlim(0.0, 90.)
ylim(-0.05, 1.0)
xticks(0:20:91.)
yticks(0:0.2:1.0)

legend()
gcf()

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

##Defining all reference values for the Medium case

RAMI_SZA = [27.,60.,83.]

RAMI_fabsRed_150_BLK =  [0.28137804, 0.46514268999999997, 0.89063486]
RAMI_frefRed_150_BLK =  [0.00923676, 0.01379672, 0.02970703]
RAMI_ftranRed_150_BLK =  [0.7093851999999999, 0.52106059, 0.07965811]

RAMI_fabsRed_150_MED =  [0.31403827, 0.49003033, 0.89432051]
RAMI_frefRed_150_MED =  [0.06195053, 0.05151941, 0.03561715]
RAMI_ftranRed_150_MED =  [0.7104761399999999, 0.52197456, 0.07977039000000001]

RAMI_fabsRed_150_SNW =  [0.5431621799999999, 0.66519762, 0.9201217300000001]
RAMI_frefRed_150_SNW =  [0.43100610000000006, 0.31581022999999997, 0.07698033]
RAMI_ftranRed_150_SNW =  [0.71754777, 0.52755972, 0.08049832999999999]
#----------------------------------------------------------------------------

##Medium case with black soil
reflRed_SZA,absRed_SZA,transRed_SZA,reflRed_clump_SZA,absRed_clump_SZA,transRed_clump_SZA=RAMI_case(1.5017, 0.0, 0.405417644);

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

SZA=0:1:85

figure(figsize=(10,5))

subplot(1,2,1)
plot(SZA, reflRed_SZA,label=["reflectance"])
plot(SZA, absRed_SZA ,label=["absorptance"])
plot(SZA, transRed_SZA,label=["transmittance"])
scatter(RAMI_SZA, RAMI_frefRed_150_BLK)
scatter(RAMI_SZA, RAMI_fabsRed_150_BLK)
scatter(RAMI_SZA, RAMI_ftranRed_150_BLK)
title("150_BLK - Default")
ylabel("Radiation Partitioning")
xlabel("Sun Zenith Angle")
xlim(0.0, 90.)
ylim(-0.05, 1.0)
xticks(0:20:91.)
yticks(0:0.2:1.0)

subplot(1,2,2)
plot(SZA, reflRed_clump_SZA,label=["reflectance"])
plot(SZA, absRed_clump_SZA ,label=["absorptance"])
plot(SZA,transRed_clump_SZA,label=["transmittance"])
scatter(RAMI_SZA, RAMI_frefRed_150_BLK,label="RAMI reflectance")
scatter(RAMI_SZA, RAMI_fabsRed_150_BLK,label="RAMI absorptance")
scatter(RAMI_SZA, RAMI_ftranRed_150_BLK,label="RAMI transmittance")
title("Clumping")
ylabel("Radiation Partitioning")
xlabel("Sun Zenith Angle")
xlim(0.0, 90.)
ylim(-0.05, 1.0)
xticks(0:20:91.)
yticks(0:0.2:1.0)

legend()
gcf()

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

##Medium case with medium soil
reflRed_SZA,absRed_SZA,transRed_SZA,reflRed_clump_SZA,absRed_clump_SZA,transRed_clump_SZA=RAMI_case(1.5017, 0.1217, 0.405417644);

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

SZA=0:1:85

figure(figsize=(10,5))

subplot(1,2,1)
plot(SZA, reflRed_SZA,label=["reflectance"])
plot(SZA, absRed_SZA ,label=["absorptance"])
plot(SZA, transRed_SZA,label=["transmittance"])
scatter(RAMI_SZA, RAMI_frefRed_150_MED)
scatter(RAMI_SZA, RAMI_fabsRed_150_MED)
scatter(RAMI_SZA, RAMI_ftranRed_150_MED)
title("150_MED - Default")
ylabel("Radiation Partitioning")
xlabel("Sun Zenith Angle")
xlim(0.0, 90.)
ylim(-0.05, 1.0)
xticks(0:20:91.)
yticks(0:0.2:1.0)

subplot(1,2,2)
plot(SZA, reflRed_clump_SZA,label=["reflectance"])
plot(SZA, absRed_clump_SZA ,label=["absorptance"])
plot(SZA,transRed_clump_SZA,label=["transmittance"])
scatter(RAMI_SZA, RAMI_frefRed_150_MED,label="RAMI reflectance")
scatter(RAMI_SZA, RAMI_fabsRed_150_MED,label="RAMI absorptance")
scatter(RAMI_SZA, RAMI_ftranRed_150_MED,label="RAMI transmittance")
title("Clumping")
ylabel("Radiation Partitioning")
xlabel("Sun Zenith Angle")
xlim(0.0, 90.)
ylim(-0.05, 1.0)
xticks(0:20:91.)
yticks(0:0.2:1.0)

legend()
gcf()

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

##Medium case with black soil
reflRed_SZA,absRed_SZA,transRed_SZA,reflRed_clump_SZA,absRed_clump_SZA,transRed_clump_SZA=RAMI_case(1.5017, 0.9640, 0.405417644);

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

SZA=0:1:85

figure(figsize=(10,5))

subplot(1,2,1)
plot(SZA, reflRed_SZA,label=["reflectance"])
plot(SZA, absRed_SZA ,label=["absorptance"])
plot(SZA, transRed_SZA,label=["transmittance"])
scatter(RAMI_SZA, RAMI_frefRed_150_SNW)
scatter(RAMI_SZA, RAMI_fabsRed_150_SNW)
scatter(RAMI_SZA, RAMI_ftranRed_150_SNW)
title("150_SNW - Default")
ylabel("Radiation Partitioning")
xlabel("Sun Zenith Angle")
xlim(0.0, 90.)
ylim(-0.05, 1.0)
xticks(0:20:91.)
yticks(0:0.2:1.0)

subplot(1,2,2)
plot(SZA, reflRed_clump_SZA,label=["reflectance"])
plot(SZA, absRed_clump_SZA ,label=["absorptance"])
plot(SZA,transRed_clump_SZA,label=["transmittance"])
scatter(RAMI_SZA, RAMI_frefRed_150_SNW,label="RAMI reflectance")
scatter(RAMI_SZA, RAMI_fabsRed_150_SNW,label="RAMI absorptance")
scatter(RAMI_SZA, RAMI_ftranRed_150_SNW,label="RAMI transmittance")
title("Clumping")
ylabel("Radiation Partitioning")
xlabel("Sun Zenith Angle")
xlim(0.0, 90.)
ylim(-0.05, 1.0)
xticks(0:20:91.)
yticks(0:0.2:1.0)

legend()
gcf()

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

##Defining all reference values for the Dense case

RAMI_SZA = [27.,60.,83.]

RAMI_fabsRed_250_BLK =  [0.46852539, 0.70426097, 0.9461774300000001]
RAMI_frefRed_250_BLK =  [0.01445858, 0.02016963, 0.03477486]
RAMI_ftranRed_250_BLK =  [0.51701603, 0.2755694, 0.01904771]

RAMI_fabsRed_250_MED =  [0.50540545, 0.72429659, 0.94742381]
RAMI_frefRed_250_MED =  [0.03953053, 0.03315039, 0.03580858]
RAMI_ftranRed_250_MED =  [0.51811911, 0.27616192, 0.01909098]

RAMI_fabsRed_250_SNW =  [0.76512258, 0.86538802, 0.9562473199999999]
RAMI_frefRed_250_SNW =  [0.21595537, 0.124503, 0.043056080000000004]
RAMI_ftranRed_250_SNW =  [0.5256125, 0.280805, 0.01935]
#----------------------------------------------------------------------------

##Dense case with black soil
reflRed_SZA,absRed_SZA,transRed_SZA,reflRed_clump_SZA,absRed_clump_SZA,transRed_clump_SZA=RAMI_case(2.5007, 0.0, 0.45946608);

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

SZA=0:1:85

figure(figsize=(10,5))

subplot(1,2,1)
plot(SZA, reflRed_SZA,label=["reflectance"])
plot(SZA, absRed_SZA ,label=["absorptance"])
plot(SZA, transRed_SZA,label=["transmittance"])
scatter(RAMI_SZA, RAMI_frefRed_250_BLK)
scatter(RAMI_SZA, RAMI_fabsRed_250_BLK)
scatter(RAMI_SZA, RAMI_ftranRed_250_BLK)
title("250_BLK - Default")
ylabel("Radiation Partitioning")
xlabel("Sun Zenith Angle")
xlim(0.0, 90.)
ylim(-0.05, 1.0)
xticks(0:20:91.)
yticks(0:0.2:1.0)

subplot(1,2,2)
plot(SZA, reflRed_clump_SZA,label=["reflectance"])
plot(SZA, absRed_clump_SZA ,label=["absorptance"])
plot(SZA,transRed_clump_SZA,label=["transmittance"])
scatter(RAMI_SZA, RAMI_frefRed_250_BLK,label="RAMI reflectance")
scatter(RAMI_SZA, RAMI_fabsRed_250_BLK,label="RAMI absorptance")
scatter(RAMI_SZA, RAMI_ftranRed_250_BLK,label="RAMI transmittance")
title("Clumping")
ylabel("Radiation Partitioning")
xlabel("Sun Zenith Angle")
xlim(0.0, 90.)
ylim(-0.05, 1.0)
xticks(0:20:91.)
yticks(0:0.2:1.0)

legend()
gcf()

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

##Medium case with medium soil
reflRed_SZA,absRed_SZA,transRed_SZA,reflRed_clump_SZA,absRed_clump_SZA,transRed_clump_SZA=RAMI_case(2.5007, 0.1217, 0.45946608);

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

SZA=0:1:85

figure(figsize=(10,5))

subplot(1,2,1)
plot(SZA, reflRed_SZA,label=["reflectance"])
plot(SZA, absRed_SZA ,label=["absorptance"])
plot(SZA, transRed_SZA,label=["transmittance"])
scatter(RAMI_SZA, RAMI_frefRed_250_MED)
scatter(RAMI_SZA, RAMI_fabsRed_250_MED)
scatter(RAMI_SZA, RAMI_ftranRed_250_MED)
title("250_MED - Default")
ylabel("Radiation Partitioning")
xlabel("Sun Zenith Angle")
xlim(0.0, 90.)
ylim(-0.05, 1.0)
xticks(0:20:91.)
yticks(0:0.2:1.0)

subplot(1,2,2)
plot(SZA, reflRed_clump_SZA,label=["reflectance"])
plot(SZA, absRed_clump_SZA ,label=["absorptance"])
plot(SZA,transRed_clump_SZA,label=["transmittance"])
scatter(RAMI_SZA, RAMI_frefRed_250_MED,label="RAMI reflectance")
scatter(RAMI_SZA, RAMI_fabsRed_250_MED,label="RAMI absorptance")
scatter(RAMI_SZA, RAMI_ftranRed_250_MED,label="RAMI transmittance")
title("Clumping")
ylabel("Radiation Partitioning")
xlabel("Sun Zenith Angle")
xlim(0.0, 90.)
ylim(-0.05, 1.0)
xticks(0:20:91.)
yticks(0:0.2:1.0)

legend()
gcf()

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

##Dense case with snowy soil
reflRed_SZA,absRed_SZA,transRed_SZA,reflRed_clump_SZA,absRed_clump_SZA,transRed_clump_SZA=RAMI_case(2.5007, 0.9640, 0.45946608);

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

SZA=0:1:85

figure(figsize=(10,5))

subplot(1,2,1)
plot(SZA, reflRed_SZA,label=["reflectance"])
plot(SZA, absRed_SZA ,label=["absorptance"])
plot(SZA, transRed_SZA,label=["transmittance"])
scatter(RAMI_SZA, RAMI_frefRed_250_SNW)
scatter(RAMI_SZA, RAMI_fabsRed_250_SNW)
scatter(RAMI_SZA, RAMI_ftranRed_250_SNW)
title("250_SNW - Default")
ylabel("Radiation Partitioning")
xlabel("Sun Zenith Angle")
xlim(0.0, 90.)
ylim(-0.05, 1.0)
xticks(0:20:91.)
yticks(0:0.2:1.0)

subplot(1,2,2)
plot(SZA, reflRed_clump_SZA,label=["reflectance"])
plot(SZA, absRed_clump_SZA ,label=["absorptance"])
plot(SZA,transRed_clump_SZA,label=["transmittance"])
scatter(RAMI_SZA, RAMI_frefRed_250_SNW,label="RAMI reflectance")
scatter(RAMI_SZA, RAMI_fabsRed_250_SNW,label="RAMI absorptance")
scatter(RAMI_SZA, RAMI_ftranRed_250_SNW,label="RAMI transmittance")
title("Clumping")
ylabel("Radiation Partitioning")
xlabel("Sun Zenith Angle")
xlim(0.0, 90.)
ylim(-0.05, 1.0)
xticks(0:20:91.)
yticks(0:0.2:1.0)

legend()
gcf()

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

# # ### Test a VZA dependence in the principal plane with clumping

## Define a few wavelengths:
wl_blue = 450.0;
wl_red = 600.0;
wl_FarRed = 740.0;
wl_Red = 685.0;
ind_wle_blue  = argmin(abs.(wl_set.wle .-wl_blue));
ind_wle_red = argmin(abs.(wl_set.wle .-wl_red));
ind_wlf_FR  = argmin(abs.(wl_set.wlf .-wl_FarRed));
ind_wlf_R  = argmin(abs.(wl_set.wlf .-wl_Red));
ind_red = argmin(abs.(wl_set.wl .-wl_Red));
ind_NIR = argmin(abs.(wl_set.wl .-800));

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

SIF_FR = Float32[]
SIF_R = Float32[]
reflVIS = Float32[]
reflNIR = Float32[]

#### Just running the code over all geometries:
##MED
soil.albedo_SW[:] .=0.1217;
#### Set sun SZA to 27 degrees
angles.tts=27.
#### Set 0 azimuth (principal plane)
angles.psi=0

##Adding clumping
canopy_rt.Ω = 1.0
#### LAI of 3:
canopy_rt.LAI = 2.5007
#### Define VZA
VZA=collect(-89.5:0.5:89.5)

for VZA_ in VZA
    angles.tto=VZA_
    compute_canopy_geometry!(canopy_rt, angles, canOpt_rt)
    compute_canopy_matrices!(arrayOfLeaves, canOpt_rt);
    simulate_short_wave!(canopy_rt, canOpt_rt, canRad_rt, sunRad_rt, soil);
    computeSIF_Fluxes!(arrayOfLeaves, canOpt_rt, canRad_rt, canopy_rt, soil, wl_set);
    #### Handpicked indices in
    push!(reflVIS, canRad_rt.alb_obs[ind_red])
    push!(reflNIR, canRad_rt.alb_obs[ind_NIR])
    push!(SIF_R , canRad_rt.SIF_obs[ind_wlf_R])
    push!(SIF_FR, canRad_rt.SIF_obs[ind_wlf_FR ])
end

##Adding clumping
canopy_rt.Ω = 0.45946608

SIF_FR_clump = Float32[]
SIF_R_clump = Float32[]
reflVIS_clump = Float32[]
reflNIR_clump = Float32[]

for VZA_ in VZA
    angles.tto=VZA_
    compute_canopy_geometry!(canopy_rt, angles, canOpt_rt)
    compute_canopy_matrices!(arrayOfLeaves, canOpt_rt);
    simulate_short_wave!(canopy_rt, canOpt_rt, canRad_rt, sunRad_rt, soil);
    computeSIF_Fluxes!(arrayOfLeaves, canOpt_rt, canRad_rt, canopy_rt, soil, wl_set);
    #### Handpicked indices in
    push!(reflVIS_clump, canRad_rt.alb_obs[ind_red])
    push!(reflNIR_clump, canRad_rt.alb_obs[ind_NIR])
    push!(SIF_R_clump , canRad_rt.SIF_obs[ind_wlf_R])
    push!(SIF_FR_clump, canRad_rt.SIF_obs[ind_wlf_FR ])
end
#----------------------------------------------------------------------------

#### Plots Visible
figure()
plot(VZA, reflVIS, color=:black,label="Red Reflectance", lw=2)
plot(VZA, SIF_R/30, color=:orange,label="Red SIF (/30)", lw=2)
plot(VZA, reflVIS_clump, color=:black, ls="--", lw=2,label="Red Reflectance w/ Clumping")
plot(VZA, SIF_R_clump/30, color=:orange, ls="--", lw=2,label="Red SIF (/30) w/ Clumping")
xlabel("Viewing Zenith Angle")
legend()
gcf()
#----------------------------------------------------------------------------

#### Plots Visible
figure()
plot(VZA, reflNIR, color=:black,label="NIR Reflectance", lw=2)
plot(VZA, SIF_FR/6, color=:orange,label="Far Red SIF (/6)", lw=2)
plot(VZA, reflNIR_clump, color=:black, ls="--", lw=2,label="NIR Reflectance w/ Clumping")
plot(VZA, SIF_FR_clump/6, color=:orange, ls="--", lw=2,label="Far Red SIF (/6) w/ Clumping")
xlabel("Viewing Zenith Angle")
legend()
gcf()
#----------------------------------------------------------------------------

# # ## BRDF sampling
#
# By going through viewing and azimuth angles, we can construct a full BRDF for reflectance and SIF emissions at different wavelengths:
#

reflVIS = Float32[]
reflNIR = Float32[]
SIF_FR = Float32[]
SIF_R  = Float32[]

##MED
soil.albedo_SW[:] .=0.1217;
angles.tts=27.
angles.psi=0
canopy_rt.LAI=2.5007
canopy_rt.Ω = 1.0
for psi=0:360
    angles.psi=psi
    for VZA=0:1:85
        angles.tto=VZA

        compute_canopy_geometry!(canopy_rt, angles, canOpt_rt)
        compute_canopy_matrices!(arrayOfLeaves, canOpt_rt);
        simulate_short_wave!(canopy_rt, canOpt_rt, canRad_rt, sunRad_rt, soil);
        computeSIF_Fluxes!(arrayOfLeaves, canOpt_rt, canRad_rt, canopy_rt, soil, wl_set);
        push!(reflVIS, canRad_rt.alb_obs[28])
        push!(reflNIR, canRad_rt.alb_obs[52])
        push!(SIF_R , canRad_rt.SIF_obs[8])
        push!(SIF_FR, canRad_rt.SIF_obs[20])
    end
end

##Adding clumping
canopy_rt.Ω = 0.45946608

SIF_FR_clump = Float32[]
SIF_R_clump = Float32[]
reflVIS_clump = Float32[]
reflNIR_clump = Float32[]

for psi=0:360
    angles.psi=psi
    for VZA=0:1:85
        angles.tto=VZA

        compute_canopy_geometry!(canopy_rt, angles, canOpt_rt)
        compute_canopy_matrices!(arrayOfLeaves, canOpt_rt);
        simulate_short_wave!(canopy_rt, canOpt_rt, canRad_rt, sunRad_rt, soil);
        computeSIF_Fluxes!(arrayOfLeaves, canOpt_rt, canRad_rt, canopy_rt, soil, wl_set);
        push!(reflVIS_clump, canRad_rt.alb_obs[28])
        push!(reflNIR_clump, canRad_rt.alb_obs[52])
        push!(SIF_R_clump , canRad_rt.SIF_obs[8])
        push!(SIF_FR_clump, canRad_rt.SIF_obs[20])
    end
end

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

A = reshape(reflNIR, ( 86,361));
B = reshape(reflVIS, ( 86,361));
SIFFER = reshape(SIF_R, ( 86,361));
SIFFER_FR = reshape(SIF_FR, ( 86,361));

A_clump = reshape(reflNIR_clump, ( 86,361));
B_clump = reshape(reflVIS_clump, ( 86,361));
SIFFER_clump = reshape(SIF_R_clump, ( 86,361));
SIFFER_FR_clump = reshape(SIF_FR_clump, ( 86,361));
#----------------------------------------------------------------------------

figure(figsize=(10,5))
subplot(1,2,1, polar=true)
grid(false)
hm = contourf(deg2rad.(collect((0:360))),collect(0:1:85),  A,  cmap=:viridis, levels=collect(0.25:0.012:0.5))
title("NIR reflectance BRDF")
yticks([])
colorbar()
subplot(1,2,2, polar=true)
grid(false)
hm = contourf(deg2rad.(collect((0:360))),collect(0:1:85),  A_clump,  cmap=:viridis, levels=collect(0.25:0.012:0.5))
title("Clumping")
yticks([])
colorbar()
gcf()

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

figure(figsize=(10,5))
subplot(1,2,1, polar=true)
grid(false)
hm = contourf(deg2rad.(collect((0:360))),collect(0:1:85),  B,  cmap=:viridis, levels=collect(0.0:0.005:0.045))
title("Red reflectance BRDF")
yticks([])
colorbar()
subplot(1,2,2, polar=true)
grid(false)
hm = contourf(deg2rad.(collect((0:360))),collect(0:1:85),  B_clump,  cmap=:viridis, levels=collect(0.0:0.005:0.045))
title("Clumping")
yticks([])
colorbar()
gcf()

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

figure(figsize=(10,5))
subplot(1,2,1, polar=true)
grid(false)
hm = contourf(deg2rad.(collect((0:360))),collect(0:1:85),  SIFFER,  cmap=:viridis, levels=collect(0.3:0.05:0.8))
title("Red SIF emission BRDF")
yticks([])
colorbar()
subplot(1,2,2, polar=true)
grid(false)
hm = contourf(deg2rad.(collect((0:360))),collect(0:1:85),  SIFFER_clump,  cmap=:viridis, levels=collect(0.3:0.05:0.8))
title("Clumping")
yticks([])
colorbar()
gcf()

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

figure(figsize=(10,5))
subplot(1,2,1, polar=true)
grid(false)
hm = contourf(deg2rad.(collect((0:360))),collect(0:1:85),  SIFFER_FR,  cmap=:viridis, levels=collect(1.:0.1:3.5))
title("Far-Red SIF emission BRDF")
yticks([])
colorbar()
subplot(1,2,2, polar=true)
grid(false)
hm = contourf(deg2rad.(collect((0:360))),collect(0:1:85),  SIFFER_FR_clump,  cmap=:viridis, levels=collect(1.:0.1:3.5))
title("Clumping")
yticks([])
colorbar()
gcf()

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

max_stars_repo_name: bbuman/Land

max_issues_repo_name: bbuman/Land

max_forks_repo_name: bbuman/Land


text:
"""
    Differentiable accelerated Ray Tracing

This is a differentiable general relativistic ray-tracer. I think that's it
"""

struct GaussianRing{R,W}
    rpeak::R
    width::W
end

@inline profile(p::GaussianRing, r) = exp(-4*log(2)*abs2( (r-p.rpeak)/p.width ))

"""
    DblPower
A double power law for a profile
"""
struct DblPower{R,A,B}
    r0::R
    a::A
    b::B
end

@inline profile(p::DblPower, r) = (r/p.r0)^p.a/(1 + (r/p.r0)^(p.a+p.b))

struct Observer{D,O}
    distance::D
    inclination::O
end

struct BAM{S, FD, V, B}
    nmax::Int
    """spectral index"""
    α::S
    """Function that computes the profile as radius `r` and redshift `z`"""
    profile::FD
    """ FluidVelocity struct"""
    β::V
    """ MagneticField struct"""
    b::B
end

function bam(nmax, spin, α, rpeak, width, βv, χ, ι, η=χ+π)
    g = Kerr(spin)
    v = FluidVelocity(βv, χ)
    b = MagneticField(ι, η)
    profile = GaussianRing(rpeak, width)
    return g, BAM{typeof(α), typeof(profile), typeof(v), typeof(b)}(nmax, α, profile, v, b)
end

struct FluidVelocity{T}
    β::T
    χ::T
    βr::T
    βϕ::T
    function FluidVelocity(β, χ)
        T = promote_type(eltype(β), eltype(χ))
        βT, χT = promote(β, χ)
        s,c = sincos(χT)
        return new{T}(βT, χT, βT*c, βT*s)
    end
end

struct MagneticField{T}
    ι::T
    η::T
    beq::T
    bz::T
    br::T
    bϕ::T
    function MagneticField(ι, η)
        T = promote_type(typeof(ι), typeof(η))
        ιT, ηT = promote(ι, η)
        beq,bz = sincos(ι)
        s,c = sincos(η)
        br = beq*c
        bϕ = beq*s
        return new{T}(ιT, ηT, beq, bz, br, bϕ)
    end
end

@inline function MagneticField(ι, v::FluidVelocity)
    return MagneticField(ι, v.χ + π)
end

@inline function magnetic_vector(b::MagneticField)
    return SVector(b.br, b.bz, b.bϕ)
end

@inline function profile(b::BAM, r)
    profile(b.profile, r)
end

@concrete struct RayCache
    interiorflag::Bool
    sβ
    λ
    η
    u₋a²
    roots
    rootdiffs
    F₀
    K₀
    k
    Ir_total
    # This is the specific stuff for interior rays
    rp
    Agl
    Bgl
    I3r
    # This is specific stuff for exterior rays
    f₀
    ffac
    Ir_turn
end

@inline function _I3(Agl, Bgl, k)
    I3r_angle = acos((Agl-Bgl)/(Agl+Bgl))
    I3r = F(I3r_angle, k)/sqrt(Agl*Bgl)
    return I3r
end

@inline function _I3total(I3r, Agl, Bgl, k, rp, r1::Real, r2::Real)
    I3rp_angle = acos( (Agl*(rp - r1) - Bgl*(rp - r2))/(Agl*(rp - r1) + Bgl*(rp - r2)) )
    I3rp = F(I3rp_angle, k)/sqrt(Agl*Bgl)
    return I3r - I3rp
end

 function interior_ray_cache(sβ, λ, η, u₋a², roots, rootdiffs, F0, K0, a)
    r1,r2,_,_ = roots
    r31, r32, r42, r41 = rootdiffs
    Agl = real(sqrt(r32*r42))
    Bgl = real(sqrt(r31*r41))
    rr1 = real(r1)
    rr2 = real(r2)
    k = ((Agl+Bgl)^2 - (rr2-rr1)^2)/(4*Agl*Bgl)
    I3r = _I3(Agl, Bgl, k)
    rp = 1+sqrt(1-a^2)
    I3total = _I3total(I3r, Agl, Bgl, k, rp, rr1, rr2)

    return RayCache(true, sβ, λ, η, u₋a², roots, rootdiffs, F0, K0,
                            k, I3total,
                            rp, Agl, Bgl, I3r,
                            zero(typeof(λ)), zero(typeof(λ)), zero(typeof(λ)) #zero this out
                    )

end

 function exterior_ray_cache(sβ, λ, η, u₋a², roots, rootdiffs, F0, K0, a)
    r31, r32, r42, r41 = rootdiffs
    k = real(r32*r41/(r31*r42))
    f₀ = F(real(asin(sqrt(r31/r41))), k)
    rp = 1+sqrt(1-a^2)
    Ir_turn = real(2/sqrt(r31*r42)*f₀)
    Ir_total = 2*Ir_turn
    ffac = 1/2*sqrt(real(r31*r42))
    T = typeof(λ)
    return RayCache(false, sβ, λ, η, u₋a², roots, rootdiffs, F0, K0,
                            k, Ir_total,
                            rp, zero(T), zero(T), zero(T),
                            f₀, ffac, Ir_turn)

end

 function raycache(sβ, λ, η, g::Kerr, o::Observer)
    a = g.spin
    up, um = get_up_um(λ, η, a)
    urat = up/um
    u₋a² = um*a^2

    roots = radialroots(λ, η, a)
    r1,r2,r3,r4 = roots
    rootdiffs = (r3-r1, r3-r2, r4-r2, r4-r1)

    F₀_sin = clamp(cos(o.inclination)/sqrt(up), -1.0, 1.0)
    F₀ = F(asin(F₀_sin), urat)
    K₀ = K(urat)
    # case 3 interior rays
    if abs(imag(r3)) > 1e-12
        return interior_ray_cache(sβ, λ, η, u₋a², roots, rootdiffs, F₀, K₀, a)
    else
        return exterior_ray_cache(sβ, λ, η, u₋a², roots, rootdiffs, F₀, K₀, a)
    end

end

function traceimg!(polim, alpha, beta, g, o, bam)
    stokesi = polim.I
    stokesq = polim.Q
    stokesu = polim.U
    @batch for C in CartesianIndices(stokesi)
        iy,ix = Tuple(C)
        x = alpha[ix]
        y = beta[iy]
        i,q,u = raytrace(x, y, g, o, bam)
        stokesi[C] = i
        stokesq[C] = q
        stokesu[C] = u
    end
end

function tracepolmap(alpha, beta, g, o, bam)
    nx,ny = length(alpha), length(beta)
    polim = StructArray((I=zeros(ny, nx), Q = zeros(ny, nx), U = zeros(ny, nx)))
    traceimg!(polim, alpha, beta, g, o, bam)
    return polim
end

function raytrace(α, β, g::Kerr, o::Observer, bam::BAM)

    # Find the geodesic charges
    λ, η = get_λ_η(α, β, o.inclination, g.spin)
    # Kill the vortical geodesics
    η < 0 && return StokesVector{typeof(λ)}(0,0,0,0)

    # precompute some stuff that will be constant over ray
    cache = raycache(sign(β), λ, η, g, o)
    # create the stokes i,q,u allocators
    T = eltype(λ)
    stokesi = zero(T); stokesq = zero(T); stokesu = zero(T)
    for n in 0:bam.nmax
        _, i, q, u = trace_nring(n, α, β, cache, g, o, bam)
        stokesi += i
        stokesq += q
        stokesu += u
    end

    return StokesVector(stokesi, stokesq, stokesu, zero(typeof(stokesi)))
end

function trace_nring(n::Int, α, β, cache::RayCache, g::Kerr, o::Observer, bam::BAM)
    r, spr = _emission_radius(n, cache)
    T = eltype(r)
    # bail out
    cache.k > 1 && return zero(T), zero(T), zero(T), zero(T)
    r < cache.rp && return zero(T), zero(T), zero(T), zero(T)
    r < zero(T) && return zero(T), zero(T), zero(T), zero(T)
    #r > 100 && return zero(T), zero(T), zero(T), zero(T)
    return _emission(n, α, β, cache.λ, cache.η, r, spr, g, o, bam)
end

function mpow(m)
    m == 0 && return one(eltype(m))
    res = -1
    for _ in 1:m
        res *= -1
    end
    return res
end

function momentum_1form(m, sβ, λ, η, spr, Δs, ℛs)
    pt = -1
    pr = clamp(spr*sqrt(abs(ℛs))/Δs, -15, 15)
    pϕ = λ
    pθ = mpow(m)*sβ*sqrt(η)
    return SVector(pt, pr, pθ, pϕ)
end

function zamo_tetrad(gs::MetricTensor{M,C}) where {M<:Kerr, C}
    Δs, Σs, Ξs, _, ωs = gs.cache
    _,r,_,_ = gs.p
    rtΞ = sqrt(Ξs)
    rtΔ = sqrt(Δs)
    ett, etϕ = 1/r*rtΞ/rtΔ, 1/r*rtΞ/rtΔ*ωs
    err = 1/r*rtΔ
    eϕϕ = r/rtΞ
    eθθ = -1/r
    T = zero(Δs)
    return @SMatrix [ ett     zero(T) zero(T) etϕ    ;
                      zero(T) err     zero(T) zero(T);
                      zero(T) zero(T) eθθ     zero(T);
                      zero(T) zero(T) zero(T) eϕϕ
                    ]
end

function penrose_walker(pcon, k, spin, rs)
    pt = pcon[1]; pr = pcon[2]; pθ = pcon[3]; pϕ = pcon[4]
    kt = k[1];    kr = k[2];    kθ = k[3];    kϕ = k[4]

    AA = (pt*kr - pr*kt) + spin*(pr*kϕ - pϕ*kr)
    BB = (rs^2 + spin^2) * (pϕ*kθ - pθ*kϕ) - spin * (pt*kθ - pθ*kt)
    κ1 = rs * AA
    κ2 = -rs * BB

    return κ1, κ2

end

# function momentum_vec(m, sβ, λ, η, spr, Δs, ℛs, r, a)
#     pt = -1/r^2*(-a*(a-λ) + (r^2 + a^2)*(r^2 + a^2 - a*λ)/Δs)
#     pr = spr*1/r^2*sqrt(ℛs)
#     pϕ = 1/r^2*(-(a-λ) + a/Δs*(r^2 + a^2 - a*λ))
#     pθ = (-1)^m*sβ*sqrt(η)/r^2
#     return SVector(pt, pr, pθ, pϕ)
# end

function _emission(n, α, β, λ, η, r, spr, g, o, bam)
    sβ = sign(β)
    m = get_nturns(n, sβ)
    #println(r, ", ", λ, ", ", η, ", ", α, ", ", β)
    # Get the metric stuff
    gs = metric(g, 0, r, π/2, 0)
    Δs, _, _, _, _ = gs.cache
    ℛs = ℛ(gs, λ, η)
    # For a 4x4
    gij = components(gs)
    # Don't worry this uses a special inverse for 4x4 matrices that should be stable
    invgij = inv(gij)
    ηij = components(metric(Minkowski(), 0, 0, 0, 0))
    # # Construct the photon momentum 1-form and vector
    pform = momentum_1form(m, sβ, λ, η, spr, Δs, ℛs)
    pcon = invgij*pform

    # # Construct the ZAMO tetrad
    ezamo = zamo_tetrad(gs)

    # # Transform to fluid frame
     vfield = bam.β
     #Λ = get_Λ(vfield.β, vfield.χ)
     #efluid = ηij*Λ*ηij*ezamo
     # -β is inverse and vectors are covariant so apply inverse transform
     Λ = get_Λ(-vfield.β, vfield.χ)
     efluid = Λ*ezamo
     pfluid = ηij*efluid*pform

    @inbounds begin
        z = 1/pfluid[1]
        lp = abs(pfluid[1]/pfluid[3])
        # #Now lets get the polarization vector
        p3 = SVector(pfluid[2], pfluid[3], pfluid[4])
    end
    b = bam.b
    bvec = magnetic_vector(b)
    f3 = cross(p3, bvec)
    f = @inbounds SVector(zero(eltype(f3)), f3[1], f3[2], f3[3])

    #Now move back to coordinate basis
    fkerr = efluid'*f
    # Constuct PW constant
    κ1, κ2 = penrose_walker(pcon, fkerr, g.spin, r)

    # screen appearance
    ν = -(α + g.spin * sin(o.inclination))
    enorm = (ν^2 + β^2)*sqrt(κ1^2+κ2^2)
    eα = (β*κ2 - ν*κ1) / enorm
    eβ = (β*κ1 + ν*κ2) / enorm
    # Get the profile value at the emission radius
    prof = profile(bam, r)*z^(3+bam.α)
    # We add a small perturbation to q and u. This prevent taking a
    # derivative of hypot at (0,0) which gives  a NaN. This is silly because this
    # occurs when there is no emission so the derivative should be zero there.
    q = -(eα^2 - eβ^2)*lp*prof + eps()
    u = -2*eα*eβ*lp*prof + eps()
    i = hypot(q, u)
    return z,i,q,u
end

@inline function get_nturns(n, sβ)
    if sβ > 0
        return sβ+n
    else
        return sβ+n+1
    end
end

@inline function minotime(m, cache)
    _Ir(cache.u₋a², m, cache.K₀, cache.sβ, cache.F₀)
end

@inline function _Ir(u₋a², m, K₀, sβ, F₀)
    return (2*m*K₀ - sβ*F₀)/sqrt(-u₋a²)
end

function _emission_radius(n::Int, cache::RayCache)
    if cache.interiorflag
        return _interior_emission_radius(n, cache)
    else
        return _exterior_emission_radius(n, cache)
    end
end

function _exterior_emission_radius(n::Int, cache)
    # get the correct number of turning point
    m = get_nturns(n, cache.sβ)
    _,_,r3,r4 = cache.roots
    r31,_,_,r41 = cache.rootdiffs

    Ir = minotime(m, cache)
    Ir > cache.Ir_total && return -one(typeof(Ir)), one(eltype(Ir))
    spr = sign(cache.Ir_turn - Ir)
    sn = Jacobi.sn(cache.ffac*Ir - cache.f₀, cache.k)
    snsqr = sn^2
    #abs(snsqr) > 1.1 && return NaN*one(eltype(real(sn))), one(eltype(real(sn)))
    r = real((r4*r31 - r3*r41*snsqr)/(r31-r41*snsqr))
    return r, spr
end

@fastmath function _interior_emission_radius(n::Int, cache)
    # get the correct number of turning point
    m = get_nturns(n, cache.sβ)
    Ir = minotime(m, cache)
    Ir > cache.Ir_total && return -one(typeof(Ir)), one(eltype(Ir))
    r1,r2,_,_ = cache.roots
    rr1 = real(r1)
    rr2 = real(r2)
    Agl = cache.Agl; Bgl = cache.Bgl
    X = sqrt(Agl*Bgl)*(Ir - cache.I3r)
    cn = Jacobi.cn(X, cache.k)
    return ((Bgl*rr2 - Agl*rr1) + (Bgl*rr2 + Agl*rr1)*cn) / ((Bgl-Agl)+(Bgl+Agl)*cn), one(typeof(X))
end