text:
##############################################################################
##
## Allocations for AllocatedLevenbergMarquardt
##
##############################################################################

type AllocatedLevenbergMarquardt{Tx1, Tx2, Ty1, Ty2} <: AbstractAllocatedOptimizer
    δx::Tx1
    dtd::Tx2
    ftrial::Ty1
    fpredict::Ty2

    function AllocatedLevenbergMarquardt(δx, dtd, ftrial, fpredict)
        length(δx) == length(dtd) || throw(DimensionMismatch("The lengths of δx and dtd must match."))
        length(ftrial) == length(fpredict) || throw(DimensionMismatch("The lengths of ftrial and fpredict must match."))
        new(δx, dtd, ftrial, fpredict)
    end
end

function AllocatedLevenbergMarquardt{Tx1, Tx2, Ty1, Ty2}(δx::Tx1, dtd::Tx2, ftrial::Ty1, fpredict::Ty2)
    AllocatedLevenbergMarquardt{Tx1, Tx2, Ty1, Ty2}(δx, dtd, ftrial, fpredict)
end

function AbstractAllocatedOptimizer{Tx, Ty, Tf, TJ, Tg}(nls::LeastSquaresProblem{Tx, Ty, Tf, TJ, Tg}, optimizer::LevenbergMarquardt)
   AllocatedLevenbergMarquardt(_zeros(nls.x), _zeros(nls.x), _zeros(nls.y), _zeros(nls.y))
end

##############################################################################
##
## Optimizer for AllocatedLevenbergMarquardt
##
##############################################################################
##############################################################################
macro levenbergtrace()
    quote
        if tracing
            update!(tr,
                    iter,
                    ssr,
                    maxabs_gr,
                    store_trace,
                    show_trace,
                    show_every
                    )
        end
    end
end

const MAX_Δ = 1e16 # minimum trust region radius
const MIN_Δ = 1e-16 # maximum trust region radius
const MIN_STEP_QUALITY = 1e-3
const GOOD_STEP_QUALITY = 0.75
const MIN_DIAGONAL = 1e-6
const MAX_DIAGONAL = 1e32

function optimize!{Tx, Ty, Tf, TJ, Tg, Toptimizer <: AllocatedLevenbergMarquardt, Tsolver}(
    anls::LeastSquaresProblemAllocated{Tx, Ty, Tf, TJ, Tg, Toptimizer, Tsolver};
            xtol::Number = 1e-8, ftol::Number = 1e-8, grtol::Number = 1e-8,
            iterations::Integer = 1_000, Δ::Number = 10.0, store_trace = false, show_trace = false, show_every = 1)

    δx, dtd = anls.optimizer.δx, anls.optimizer.dtd
    ftrial, fpredict = anls.optimizer.ftrial, anls.optimizer.fpredict
    x, fcur, f!, J, g! = anls.x, anls.y, anls.f!, anls.J, anls.g!

    decrease_factor = 2.0
    # initialize
    f_calls,  g_calls, mul_calls = 0, 0, 0
    converged, x_converged, f_converged, gr_converged, converged =
        false, false, false, false, false
    f!(x, fcur)
    f_calls += 1
    ssr = sumabs2(fcur)
    maxabs_gr = Inf
    need_jacobian = true

    eTx, eTy = eltype(x), eltype(fcur)

    iter = 0

    tr = OptimizationTrace()
    tracing = store_trace || show_trace
    @levenbergtrace

    while !converged && iter < iterations
        iter += 1
        check_isfinite(x)

        # compute step
        if need_jacobian
            g!(x, J)
            g_calls += 1
            need_jacobian = false
        end
        colsumabs2!(dtd, J)
        clamp!(dtd, MIN_DIAGONAL, MAX_DIAGONAL)
        scale!(dtd, 1/Δ)
        δx, lmiter = A_ldiv_B!(δx, J, fcur, dtd,  anls.solver)
        mul_calls += lmiter
        #update x
        axpy!(-one(eTx), δx, x)
        f!(x, ftrial)
        f_calls += 1

        # trial ssr
        trial_ssr = sumabs2(ftrial)

        # predicted ssr
        A_mul_B!(one(eTx), J, δx, zero(eTx), fpredict)
        mul_calls += 1
        axpy!(-one(eTy), fcur, fpredict)
        predicted_ssr = sumabs2(fpredict)
        ρ = (ssr - trial_ssr) / (ssr - predicted_ssr)

        Ac_mul_B!(one(eTx), J, fcur, zero(eTx), dtd)
        maxabs_gr = maxabs(dtd)
        mul_calls += 1

        x_converged, f_converged, gr_converged, converged =
            assess_convergence(δx, x, maxabs_gr, ssr, trial_ssr, xtol, ftol, grtol)

        if ρ > MIN_STEP_QUALITY
            copy!(fcur, ftrial)
            ssr = trial_ssr
            # increase trust region radius (from Ceres solver)
            Δ = min(Δ / max(1/3, 1.0 - (2.0 * ρ - 1.0)^3), MAX_Δ)
            decrease_factor = 2.0
            need_jacobian = true
        else
            # revert update
            axpy!(one(eTx), δx, x)
            Δ = max(Δ / decrease_factor , MIN_Δ)
            decrease_factor *= 2.0
        end
        @levenbergtrace
    end
    LeastSquaresResult("LevenbergMarquardt", x, ssr, iter, converged,
                        x_converged, xtol, f_converged, ftol, gr_converged, grtol, tr,
                        f_calls, g_calls, mul_calls)
end


text:
@doc """
### SimpleODEMod.calc_x2_4t3

  Calculates and returns x^2 + 4*t^3
"""->
function calc_x2_4t3(coords::AbstractArray{Tmsh}, params::ParamType2, t) where Tmsh

  x = coords[1]
  y = coords[2]

  u = x^2 + 4*(t^3)

  return u
end

@doc """
### SimpleODEMod.calc_x2_3t2

  Calculates and returns x^2 + 3*t^2
"""->
function calc_x2_3t2(coords::AbstractArray{Tmsh}, params::ParamType2, t) where Tmsh

  x = coords[1]
  y = coords[2]

  u = x^2 + 3*(t^2)

  return u
end

@doc """
### SimpleODEMod.calc_x2_t2

  Calculates and returns x^2 + t^2
"""->
function calc_x2_t2(coords::AbstractArray{Tmsh}, params::ParamType2, t) where Tmsh

  x = coords[1]
  y = coords[2]

  u = x^2 + t^2

  return u
end

@doc """
### SimpleODEMod.calc_4t3

  Calculates and returns 4t^3
  (du/dt of calc_x2_t4)
"""->
function calc_4t3(coords::AbstractArray{Tmsh}, params::ParamType2, t) where Tmsh

  x = coords[1]
  y = coords[2]

  u = 4*(t^3)

  return u
end

@doc """
### SimpleODEMod.calc_3t2

  Calculates and returns 3t^2
  (du/dt of calc_x2_t3)
"""->
function calc_3t2(coords::AbstractArray{Tmsh}, params::ParamType2, t) where Tmsh

  x = coords[1]
  y = coords[2]

  u = 3*(t^2)

  return u
end

@doc """
### SimpleODEMod.calc_2t

  Calculates and returns 2t
  (du/dt of calc_x2_t2)
"""->
function calc_2t(coords::AbstractArray{Tmsh}, params::ParamType2, t) where Tmsh

  x = coords[1]
  y = coords[2]

  u = 2*t

  return u
end

ext: jl


text:
# -*- coding: utf-8 -*-
# ---
# jupyter:
#   jupytext:
#     text_representation:
#       extension: .jl
#       format_name: light
#       format_version: '1.3'
#       jupytext_version: 0.8.6
#   kernelspec:
#     display_name: Julia 1.0.3
#     language: julia
#     name: julia-1.0
# ---

module DiscreteTimeSIR
import Random
using RandomNumbers
using DataFrames

@inline @fastmath function randbn(n,p,rng)
    q = 1.0 - p
    s = p/q
    a = (n+1)*s
    r = exp(n*log(q))
    x = 0
    u = rand()
    while true
        if (u < r)
            return x
        end
        u -= r
        x += 1
        r *= (a/x)-s
    end
end

@inline @fastmath function sir(u, prams, rng)

    (S, I, R, Y) = u
    (β, γ, ι, N, δt) = prams
    λ = β * (I + ι) / N
    ifrac = 1.0 - exp(-λ * δt)
    rfrac = 1.0 - exp(-γ * δt)
    infection = randbn(S, ifrac, rng)
    recovery = randbn(I, rfrac, rng)
    return (S - infection, I + infection - recovery, R + recovery, Y + infection)
end

function simulate(r)
    prams = (0.1, 0.05, 0.01, 1000.0, 0.1)
    tf = 200
    t = 0:0.1:tf
    tl = length(t)
    S = zeros(tl)
    I = zeros(tl)
    R = zeros(tl)
    Y = zeros(tl)
    u0 = (999, 1, 0, 0)
    (S[1],I[1],R[1],Y[1]) = u0
    u = u0
    for i in 2:tl
        u = sir(u, prams, r)
        (S[i],I[i],R[i],Y[i]) = u
    end
    return DataFrame(Time=t,S=S,I=I,R=R,Y=Y)
end

# NOTE: Removed GTRI's fancy number generator
# seed = 42
# r = Xorshifts.Xorshift128Plus(seed);

Random.seed!(42);

# NOTE: printing output for potential correctness verification
sir_out = simulate(r);
print(sir_out)

end


text:
#2-dim version

include("grid.jl")
include("boundaries.jl")
include("visualize.jl")
include("eos.jl")
include("reconstruct.jl")
include("rsolvers.jl")
include("gravity.jl")
include("viscosity.jl")
include("integrators.jl")

#using Winston
using Plots

#Primary variables to conserved variables
function prim2con(rho::AbstractMatrix,
                  velx::AbstractMatrix,
                  vely::AbstractMatrix,
                  eps::AbstractMatrix)

    q = zeros(size(rho, 1), size(rho, 2), 4)
    q[:, :, 1] = rho
    q[:, :, 2] = rho .* velx
    q[:, :, 3] = rho .* vely
    q[:, :, 4] = rho .* eps .+ 0.5rho .* (velx.^2.0 .+ vely.^2.0)

    return q
end

#Conserved variables to primary variables
function con2prim(q)
    rho = q[:, :, 1]
    velx = q[:, :, 2] ./ rho
    vely = q[:, :, 3] ./ rho
    eps = q[:, :, 4] ./ rho .- 0.5(velx.^2.0 .+ vely.^2.0)

    press = eos_press(rho, eps, gamma)

    return rho, eps, press, velx, vely
end

#time step calculation
function calc_dt(hyd, dtp)
    cs = sqrt(eos_cs2(hyd.rho, hyd.eps, gamma))
    dtnew = 1.0
    for j = (hyd.g+1):(hyd.ny-hyd.g+1), i = (hyd.g+1):(hyd.nx-hyd.g+1)
        dtnew = min(dtnew, (hyd.x[i+1] - hyd.x[i]) / max(abs(hyd.velx[j, i]+cs[j, i]), abs(hyd.velx[j, i]-cs[j, i])))
        dtnew = min(dtnew, (hyd.y[j+1] - hyd.y[j]) / max(abs(hyd.vely[j, i]+cs[j, i]), abs(hyd.vely[j, i]-cs[j, i])))
    end

    dtnew = min(cfl*dtnew, 1.05*dtp)

    return dtnew
end

#snapshot of the simulation
function snapshot(hyd, var, name, it)

    #cm = Uint32[Color.convert(Color.RGB24,c) for c in flipud(Color.colormap("RdBu"))]
    xs = hyd.g+1
    xe = hyd.nx-hyd.g
    ye = hyd.ny-hyd.g

    hdata = var[xs:ye, xs:xe]
    #pf=FramedPlot()
    clims = (0., 2.5)
    #img = Winston.data2rgb(hdata, clims, cm)
	  #img = Winston.imagesc(hdata)
    #add(pf, Image((hyd.x[xs], hyd.x[xe]), (hyd.y[xs], hyd.y[ye]), img;))
		pf = heatmap(hdata,ratio=:equal)
		#add(pf, img)
    #setattr(pf, xrange=(hyd.x[xs], hyd.x[xe]))
    #setattr(pf, yrange=(hyd.y[xs], hyd.y[ye]))
    #setattr(pf, title="rho")

    #setattr(pf.frame,
    #        draw_spine=false,
    #        draw_ticks=false,
    #        draw_ticklabels=false)

    numb = lpad(it, 3, '0')
    savefig(pf,"$(name)_$(numb).png")
		closeall()

    return nothing
end

#Construct the right hand side of the group of equations
function calc_rhs(hyd, dt)

    #reconstruction and prim2con
    hyd = reconstruct(hyd)

    #apply source terms for the interfaces
    #hyd = source_terms(hyd, dt)

    #unsplit flux
    fx, fy = uflux(hyd, dt)

    #add artificial viscosity
    #vfx, vfy = artificial_viscosity(hyd, dt)
    #fx += vfx
    #fy += vfy

    #flux difference
    fluxdiff = zeros(hyd.ny, hyd.nx, 4)
    for j = (hyd.g-1):(hyd.ny-hyd.g+1), i = (hyd.g-1):(hyd.nx-hyd.g+1)
        dx = abs(hyd.x[i+1]-hyd.x[i])
        dy = abs(hyd.y[j+1]-hyd.y[j])

        for k = 1:4
            fluxdiff[j, i, k] = (fx[j, i, k] - fx[j, i-1, k])/dx + (fy[j, i, k] - fy[j-1, i, k])/dy
        end
    end

    return -fluxdiff
end

###############
# main program
###############

function evolve(hyd, tend, gamma, cfl, nx, ny)

    time = range(0.0, length=100, stop=tend)
    it = 1

    dt = 1.0e-8
    dtp = dt

    #get initial timestep
    dt = calc_dt(hyd, dt)

    #initial prim2con
    hyd.q = prim2con(hyd.rho, hyd.velx, hyd.vely, hyd.eps)

    t = 0.0
    i = 0

    visualize(hyd)

    while t < tend

        if i % 10 == 0
            println("$i $t $dt")
            #visualize(hyd)
        end

        if time[it] <= t
            snapshot(hyd, hyd.rho, "density", it)
            snapshot(hyd, hyd.press, "pressure", it)
            it += 1
        end

        #calculate new timestep
        dt = calc_dt(hyd, dt)

        #time integration
        hyd = RK2(hyd, dt)

        #update time
        t += dt
        i += 1

    end

    return hyd
end

#basic parameters
gamma = 1.4
cfl = 0.4

nx = 50
ny = 50
tend = 1.0e-1

#initialize
hyd = data2d(nx, ny)
#r32x, r32y = r32kernel(hyd.x, hyd.y)

#set up grid

include("../tests/1dim_hd.jl")
include("../tests/2dim_hd.jl")

hyd = grid_setup(hyd, 0.0, 1.0, 0.0, 1.0)
#hyd = grid_setup(hyd, 0.0, 1.0, 0.0, 0.5)
#hyd = grid_setup(hyd, -0.25, 0.25, -0.75, 0.75)
#hyd = grid_setup(hyd, 0., 1., -0.5, 0.5)

#set up initial data
#hyd = setup_taylor2(hyd)
hyd = setup_blast(hyd)
#hyd = setup_tubexy(hyd)
#hyd = setup_tubey(hyd)
#hyd = setup_tubex(hyd)
#hyd = setup_collision(hyd)
#hyd = setup_fall(hyd)
#hyd = setup_kh(hyd)

#visualize(hyd)

#main integration loop
hyd = evolve(hyd, tend, gamma, cfl, nx, ny)
#hyd = @profile evolve(hyd, tend, gamma, cfl, nx, ny)

max_stars_repo_name: gzahl/hydro

max_issues_repo_name: gzahl/hydro

max_forks_repo_name: gzahl/hydro


text:
module BioFeatures

# Region implies a range or something, Feature could be a single point

export FeatureTable
export Features
export FeaturesSlim
export add_dataframe_to_features!
export insertintervals!
export insertfeatures!
include(Pkg.dir("BioFeatures","src","Features.jl"))

end # module


text:
solaroptions() = Dict(
    :gisregion => "Europe8",            # "Europe8", "Eurasia38", "Scand3"
    :filenamesuffix => "",              # e.g. "_landx2" to save high land availability data as "GISdata_solar2018_Europe8_landx2.mat"

    :pv_density => 45,                  # Solar PV land use 45 Wp/m2 = 45 MWp/km2 (includes PV efficiency & module spacing, add latitude dependency later)
    :csp_density => 35,                 # CSP land use 35 W/m2

    :pvroof_area => .05,                # area available for rooftop PV after the masks have been applied
    :plant_area => .05,                 # area available for PV or CSP plants after the masks have been applied

    :distance_elec_access => 150,       # max distance to grid [km] (for solar classes of category B)
    :plant_persons_per_km2 => 150,      # not too crowded, max X persons/km2 (both PV and CSP plants)
    :pvroof_persons_per_km2 => 200,     # only in populated areas, so AT LEAST x persons/km2
                                        # US census bureau requires 1000 ppl/mile^2 = 386 ppl/km2 for "urban" (half in Australia)
                                        # roughly half the people of the world live at density > 300 ppl/km2
    :exclude_landtypes => [0,1,2,3,4,5,8],       # exclude water and forests. See codes in table below.
    :protected_codes => [1,2,3,4,5,8],  # IUCN codes to be excluded as protected areas. See codes in table below.

    :scenarioyear => "ssp2_2050",       # default scenario and year for population and grid access datasets
    :era_year => 2018,                  # which year of the ERA5 time series to use

    :res => 0.01,                       # resolution of auxiliary datasets [degrees per pixel]
    :erares => 0.28125,                 # resolution of ERA5 datasets [degrees per pixel]

    :pvclasses_min => [0.08,0.14,0.18,0.22,0.26],   # lower bound on annual PV capacity factor for class X    [0:0.01:0.49;]
    :pvclasses_max => [0.14,0.18,0.22,0.26,1.00],   # upper bound on annual PV capacity factor for class X    [0.01:0.01:0.50;]
    :cspclasses_min => [0.10,0.18,0.24,0.28,0.32],  # lower bound on annual CSP capacity factor for class X
    :cspclasses_max => [0.18,0.24,0.28,0.32,1.00],  # upper bound on annual CSP capacity factor for class X

    :downsample_masks => 1,     # set to 2 or higher to scale down mask sizes to avoid GPU errors in Makie plots for large regions
    :classB_threshold => 0.001  # minimum share of pixels within distance_elec_access km that must have grid access
                                # for a pixel to be considered for solar class B.
)
    # Land types
    #     0      'Water'
    #     1      'Evergreen Needleleaf Forests'
    #     2      'Evergreen Broadleaf Forests'
    #     3      'Deciduous Needleleaf Forests'
    #     4      'Deciduous Broadleaf Forests'
    #     5      'Mixed Forests'
    #     6      'Closed Shrublands'
    #     7      'Open Shrublands'
    #     8      'Woody Savannas'
    #     9      'Savannas'
    #    10      'Grasslands'
    #    11      'Permanent Wetlands'
    #    12      'Croplands'
    #    13      'Urban'
    #    14      'Cropland/Natural'
    #    15      'Snow/Ice'
    #    16      'Barren'

    # Protected areas (IUCN codes from the WDPA)
    #    1      'Ia'                'Strict Nature Reserve'
    #    2      'Ib'                'Wilderness Area'
    #    3      'II'                'National Park'
    #    4      'III'               'Natural Monument'
    #    5      'IV'                'Habitat/Species Management'
    #    6      'V'                 'Protected Landscape/Seascape'
    #    7      'VI'                'Managed Resource Protected Area'
    #    8      'Not Reported'      'Not Reported'
    #    9      'Not Applicable'    'Not Applicable'
    #    10     'Not Assigned'      'Not Assigned'

mutable struct SolarOptions
    gisregion               ::String
    filenamesuffix          ::String
    pv_density              ::Float64           # W/m2
    csp_density             ::Float64           # W/m2
    pvroof_area             ::Float64           # share [0-1]
    plant_area              ::Float64           # share [0-1]
    distance_elec_access    ::Float64           # km
    plant_persons_per_km2   ::Float64           # persons/km2
    pvroof_persons_per_km2  ::Float64           # persons/km2
    exclude_landtypes       ::Vector{Int}
    protected_codes         ::Vector{Int}
    scenarioyear            ::String
    era_year                ::Int
    res                     ::Float64           # degrees/pixel
    erares                  ::Float64           # degrees/pixel
    pvclasses_min           ::Vector{Float64}
    pvclasses_max           ::Vector{Float64}
    cspclasses_min          ::Vector{Float64}
    cspclasses_max          ::Vector{Float64}
    downsample_masks        ::Int
    classB_threshold        ::Float64
end

SolarOptions() = SolarOptions("","",0,0,0,0,0,0,0,[],[],"",0,0,0,[],[],[],[],0,0.0)

function SolarOptions(d::Dict{Symbol,Any})
    options = SolarOptions()
    for (key,val) in d
        setproperty!(options, key, val)
    end
    return options
end

function GISsolar(; savetodisk=true, plotmasks=false, optionlist...)

    # IMPORTANT!! The function makesolarera5() uses ERA5 solar datasets to
    # calculate Global Tilted Irradiance (GTI) for solar PV and Direct Normal
    # Irradiance (DNI) for CSP solar towers. If we want CSP parabolic troughs
    # then we need to add that dataset in makesolarera5() (and capacity factors
    # will become somewhat lower).

    # NOTE ON SOLAR UNIT: the solar irradiance sets are in kW/m2. Since the
    # irradiance value used to represent "standard testing conditions" for PV
    # is 1000 W/m2, the solar datasets also directly give the capacity factor.
    # Actual insolation can occasionally go above 1000 W/m2.

    # Ideally, we should make direct assumptions of PV module efficiency as a
    # function of air temperature (see e.g. Bett & Thornton appendix A2), but
    # for now efficiency is included in our assumption of :pv_density. Wind
    # speed also affects PV module temperature and efficiency. However, the
    # uncertainties in :pv_density and :plant_area are so large that efficiency
    # variations as a function of temperature don't matter.

    options = SolarOptions(merge(solaroptions(), optionlist))
    @unpack gisregion, era_year, filenamesuffix, pv_density, csp_density, downsample_masks = options

    regions, offshoreregions, regionlist, gridaccess, popdens, topo, land, protected, lonrange, latrange =
                read_datasets(options)

    mask_rooftop, mask_plantA, mask_plantB =
        create_solar_masks(options, regions, gridaccess, popdens, land, protected, lonrange, latrange,
                            plotmasks=plotmasks, downsample=downsample_masks)

    plotmasks == :onlymasks && return nothing

    meanGTI, solarGTI, meanDNI, solarDNI = read_solar_datasets(options, lonrange, latrange)

    CF_pvrooftop, CF_pvplantA, CF_pvplantB, CF_cspplantA, CF_cspplantB, solar_overlap_areaA, solar_overlap_areaB,
            capacity_pvrooftop, capacity_pvplantA, capacity_pvplantB, capacity_cspplantA, capacity_cspplantB =
        calc_solar_vars(options, meanGTI, solarGTI, meanDNI, solarDNI, regions, offshoreregions, regionlist,
                mask_rooftop, mask_plantA, mask_plantB, lonrange, latrange)

    if savetodisk
        mkpath(in_datafolder("output"))
        matopen(in_datafolder("output", "GISdata_solar$(era_year)_$gisregion$filenamesuffix.mat"), "w", compress=true) do file
            write(file, "CFtime_pvrooftop", CF_pvrooftop)
            write(file, "CFtime_pvplantA", CF_pvplantA)
            write(file, "CFtime_pvplantB", CF_pvplantB)
            write(file, "CFtime_cspplantA", CF_cspplantA)
            write(file, "CFtime_cspplantB", CF_cspplantB)
            write(file, "capacity_pvrooftop", capacity_pvrooftop)
            write(file, "capacity_pvplantA", capacity_pvplantA)
            write(file, "capacity_pvplantB", capacity_pvplantB)
            write(file, "capacity_cspplantA", capacity_cspplantA)
            write(file, "capacity_cspplantB", capacity_cspplantB)
            write(file, "solar_overlap_areaA", solar_overlap_areaA)
            write(file, "solar_overlap_areaB", solar_overlap_areaB)
            write(file, "pv_density", pv_density)
            write(file, "csp_density", csp_density)
        end
    end

    nothing
    # return CF_pvrooftop, CF_pvplantA, CF_pvplantB, CF_cspplantA, CF_cspplantB,
    #         capacity_pvrooftop, capacity_pvplantA, capacity_pvplantB, capacity_cspplantA, capacity_cspplantB
end

function read_solar_datasets(options, lonrange, latrange)
    @unpack res, erares, era_year = options

    println("Reading ERA5 solar datasets...")
    eralonranges, eralatrange = eraranges(lonrange, latrange, res, erares)

    @time meanGTI, solarGTI, meanDNI, solarDNI = h5open(in_datafolder("era5solar$era_year.h5"), "r") do file
        if length(eralonranges) == 1
            file["meanGTI"][eralonranges[1], eralatrange],
                file["GTI"][:,eralonranges[1], eralatrange],
                file["meanDNI"][eralonranges[1], eralatrange],
                file["DNI"][:,eralonranges[1], eralatrange]
        else
            [file["meanGTI"][eralonranges[1], eralatrange]; file["meanGTI"][eralonranges[2], eralatrange]],
                [file["GTI"][:, eralonranges[1], eralatrange] file["GTI"][:, eralonranges[2], eralatrange]],
                [file["meanDNI"][eralonranges[1], eralatrange]; file["meanDNI"][eralonranges[2], eralatrange]],
                [file["DNI"][:, eralonranges[1], eralatrange] file["DNI"][:, eralonranges[2], eralatrange]]
        end
    end
    return meanGTI, solarGTI, meanDNI, solarDNI
end

function create_solar_masks(options, regions, gridaccess, popdens, land, protected, lonrange, latrange; plotmasks=false, downsample=1)
    @unpack res, gisregion, exclude_landtypes, protected_codes, distance_elec_access, plant_persons_per_km2,
            pvroof_persons_per_km2, classB_threshold, filenamesuffix = options

    println("Creating masks...")

    goodland = (regions .> 0)
    for i in exclude_landtypes
        goodland[land .== i] .= false
    end
    protected_area = zeros(Bool, size(protected))
    for i in protected_codes
        protected_area[protected .== i] .= true
    end

    # Pixels with electricity access for onshore wind A
    gridA = (gridaccess .> 0)

    # Pixels with electricity access for onshore wind B and offshore wind
    km_per_degree = π*2*6371/360
    disk = diskfilterkernel(distance_elec_access/km_per_degree/res)
    gridB = (imfilter(gridaccess, disk) .> max(1e-9, classB_threshold)) # avoid artifacts if classB_threshold == 0

    # println("MAKE SURE MASKS DON'T OVERLAP! (regions & offshoreregions, mask_*)")

    # all mask conditions
    mask_rooftop = gridA .& (popdens .> pvroof_persons_per_km2) .& .!protected_area
    mask_plantA = gridA .& (popdens .< plant_persons_per_km2) .& goodland .& .!protected_area
    mask_plantB = (gridB .& .!gridA) .& (popdens .< plant_persons_per_km2) .& goodland .& .!protected_area

    if plotmasks != false   # can == :onlymasks as well
        # drawmap(land)
        isregion = (regions .> 0) .& (regions .!= NOREGION)

        # mask values refer to colors in ColorBrewer Set2_7:
        # https://juliagraphics.github.io/ColorSchemes.jl/stable/basics/#colorbrewer-1
        masks = zeros(Int16, size(regions))
        masks[(masks .== 0) .& (popdens .> plant_persons_per_km2)] .= 2
        masks[(masks .== 0) .& protected_area] .= 3
        masks[(masks .== 0) .& .!gridA .& .!gridB] .= 4
        masks[(masks .== 0) .& .!goodland] .= 1
        masks[(masks .== 0) .& .!gridA .& gridB] .= 6
        masks[(masks .== 0) .& isregion] .= 5
        masks[regions .== 0] .= 0
        masks[regions .== NOREGION] .= NOREGION
        legendtext = ["bad land type", "high population", "protected area", "no grid", "solar plant A", "solar plant B", "", ""]
        maskmap("$(gisregion)_masks_solar$filenamesuffix", masks, legendtext, lonrange, latrange; legend=true, downsample=downsample)
    end

    return mask_rooftop, mask_plantA, mask_plantB
end

function increment_solarCF!(cf::AbstractVector{<:AbstractFloat}, solardata::AbstractVector{<:AbstractFloat})
    @inbounds for i = 1:length(cf)
        cf[i] += solardata[i]
    end
end

function makesolarclasses(options, meanGTI, meanDNI)
    println("Allocating pixels to classes using ERA5 annual means...")

    @unpack pvclasses_min, pvclasses_max, cspclasses_min, cspclasses_max = options

    pvclass = getclasses(meanGTI, pvclasses_min, pvclasses_max)
    cspclass = getclasses(meanDNI, cspclasses_min, cspclasses_max)

    return pvclass, cspclass
end

function calc_solar_vars(options, meanGTI, solarGTI, meanDNI, solarDNI, regions, offshoreregions, regionlist,
                mask_rooftop, mask_plantA, mask_plantB, lonrange, latrange)

    pvclass, cspclass = makesolarclasses(options, meanGTI, meanDNI)
    eralons, eralats, lonmap, latmap, cellarea = eralonlat(options, lonrange, latrange)

    println("Calculating GW potential and hourly capacity factors for each region and class...")
    # println("Interpolate ERA5 insolation later (maybe 4x runtime).")

    @unpack era_year, pvclasses_min, cspclasses_min, res, erares, pv_density, csp_density, pvroof_area, plant_area = options

    numreg = length(regionlist)
    npvclasses, ncspclasses = length(pvclasses_min), length(cspclasses_min)
    yearlength, nlons, nlats = size(solarGTI)
    firsttime = DateTime(era_year, 1, 1)

    capacity_pvrooftop = zeros(numreg,npvclasses)
    capacity_pvplantA = zeros(numreg,npvclasses)
    capacity_pvplantB = zeros(numreg,npvclasses)
    capacity_cspplantA = zeros(numreg,ncspclasses)
    capacity_cspplantB = zeros(numreg,ncspclasses)
    CF_pvrooftop = zeros(yearlength,numreg,npvclasses)
    CF_pvplantA = zeros(yearlength,numreg,npvclasses)
    CF_pvplantB = zeros(yearlength,numreg,npvclasses)
    CF_cspplantA = zeros(yearlength,numreg,ncspclasses)
    CF_cspplantB = zeros(yearlength,numreg,ncspclasses)
    count_pvrooftop = zeros(Int,numreg,npvclasses)
    count_pvplantA = zeros(Int,numreg,npvclasses)
    count_pvplantB = zeros(Int,numreg,npvclasses)
    count_cspplantA = zeros(Int,numreg,ncspclasses)
    count_cspplantB = zeros(Int,numreg,ncspclasses)
    solar_overlap_areaA = zeros(numreg,npvclasses,ncspclasses)
    solar_overlap_areaB = zeros(numreg,npvclasses,ncspclasses)

    # Run times vary wildly depending on geographical area (because of far offshore regions with mostly zero wind speeds).
    # To improve the estimated time of completing the progress bar, iterate over latitudes in random order.
    Random.seed!(1)
    updateprogress = Progress(nlats, 1)
    for j in randperm(nlats)
        eralat = eralats[j]
        colrange = latmap[lat2col(eralat+erares/2, res):lat2col(eralat-erares/2, res)-1]
        for i = 1:nlons
            meanGTI[i,j] == 0 && meanDNI[i,j] == 0 && continue
            GTI = solarGTI[:, i, j]
            DNI = solarDNI[:, i, j]
            eralon = eralons[i]
            # get all high resolution row and column indexes within this ERA5 cell
            rowrange = lonmap[lon2row(eralon-erares/2, res):lon2row(eralon+erares/2, res)-1]

            for c in colrange, r in rowrange
                (c == 0 || r == 0) && continue
                reg = regions[r,c]
                (reg == 0 || reg == NOREGION) && continue

                area = cellarea[c]
                class = pvclass[i,j]
                # can't use elseif here, probably some overlap in the masks
                # @views is needed to make sure increment_windCF!() works with matrix slices
                # also faster since it avoids making copies
                @views if class > 0
                    if mask_rooftop[r,c] > 0
                        capacity_pvrooftop[reg,class] += 1/1000 * pv_density * pvroof_area * area
                        increment_solarCF!(CF_pvrooftop[:,reg,class], GTI)
                        count_pvrooftop[reg,class] += 1
                    elseif mask_plantA[r,c] > 0
                        capacity_pvplantA[reg,class] += 1/1000 * pv_density * plant_area * area
                        increment_solarCF!(CF_pvplantA[:,reg,class], GTI)
                        count_pvplantA[reg,class] += 1
                    elseif mask_plantB[r,c] > 0
                        capacity_pvplantB[reg,class] += 1/1000 * pv_density * 2 * plant_area * area
                        increment_solarCF!(CF_pvplantB[:,reg,class], GTI)
                        count_pvplantB[reg,class] += 1
                    end
                end

                class_pv = class
                class = cspclass[i,j]
                # @views is needed to make sure increment_windCF!() works with matrix slices
                # also faster since it avoids making copies
                @views if class > 0
                    if mask_plantA[r,c] > 0
                        capacity_cspplantA[reg,class] += 1/1000 * csp_density * plant_area * area
                        increment_solarCF!(CF_cspplantA[:,reg,class], DNI)
                        count_cspplantA[reg,class] += 1
                    elseif mask_plantB[r,c] > 0
                        capacity_cspplantB[reg,class] += 1/1000 * csp_density * 2 * plant_area * area
                        increment_solarCF!(CF_cspplantB[:,reg,class], DNI)
                        count_cspplantB[reg,class] += 1
                    end
                end

                if class_pv > 0 && class > 0
                    if mask_plantA[r,c] > 0
                        solar_overlap_areaA[reg,class_pv,class] += 1/1000 * plant_area * area
                    elseif mask_plantB[r,c] > 0
                        solar_overlap_areaB[reg,class_pv,class] += 1/1000 * 2 * plant_area * area
                    end
                end
            end
        end
        next!(updateprogress)
    end

    for y = 1:yearlength
        CF_pvrooftop[y,:,:] ./= count_pvrooftop
        CF_pvplantA[y,:,:] ./= count_pvplantA
        CF_pvplantB[y,:,:] ./= count_pvplantB
        CF_cspplantA[y,:,:] ./= count_cspplantA
        CF_cspplantB[y,:,:] ./= count_cspplantB
    end

    return CF_pvrooftop, CF_pvplantA, CF_pvplantB, CF_cspplantA, CF_cspplantB, solar_overlap_areaA, solar_overlap_areaB,
            capacity_pvrooftop, capacity_pvplantA, capacity_pvplantB, capacity_cspplantA, capacity_cspplantB
end

# Quick and ugly copy/paste hack to create resource maps for solar classes combined with masks.
function GISsolarmap(; optionlist...)
    options = SolarOptions(merge(solaroptions(), optionlist))
    @unpack gisregion, era_year, filenamesuffix = options

    regions, offshoreregions, regionlist, gridaccess, popdens, topo, land, protected, lonrange, latrange =
                read_datasets(options)
    meanGTI, solarGTI, meanDNI, solarDNI = read_solar_datasets(options, lonrange, latrange)

    mask_rooftop, mask_plantA, mask_plantB =
        create_solar_masks(options, regions, gridaccess, popdens, land, protected, lonrange, latrange, plotmasks=true)

    pvmap, pvrooftopmap, cspmap =
        calc_solar_map(options, meanGTI, solarGTI, meanDNI, solarDNI, regions, offshoreregions, regionlist,
                mask_rooftop, mask_plantA, mask_plantB, lonrange, latrange)

    return pvmap, pvrooftopmap, cspmap
end

function calc_solar_map(options, meanGTI, solarGTI, meanDNI, solarDNI, regions, offshoreregions, regionlist,
                mask_rooftop, mask_plantA, mask_plantB, lonrange, latrange)

    pvclass, cspclass = makesolarclasses(options, meanGTI, meanDNI)
    eralons, eralats, lonmap, latmap, cellarea = eralonlat(options, lonrange, latrange)

    println("Calculating GW potential and hourly capacity factors for each region and class...")
    println("Interpolate ERA5 insolation later (maybe 4x runtime).")

    @unpack era_year, pvclasses_min, cspclasses_min, res, erares, pv_density, csp_density, pvroof_area, plant_area = options

    numreg = length(regionlist)
    npvclasses, ncspclasses = length(pvclasses_min), length(cspclasses_min)
    yearlength, nlons, nlats = size(solarGTI)
    firsttime = DateTime(era_year, 1, 1)

    # pvmap = zeros(size(regions))
    # pvrooftopmap = zeros(size(regions))
    # cspmap = zeros(size(regions))
    pvmap = zeros(Int16, size(regions))
    pvrooftopmap = zeros(Int16, size(regions))
    cspmap = zeros(Int16, size(regions))

    # Run times vary wildly depending on geographical area (because of far offshore regions with mostly zero wind speeds).
    # To improve the estimated time of completing the progress bar, iterate over latitudes in random order.
    Random.seed!(1)
    updateprogress = Progress(nlats, 1)
    @inbounds for j in randperm(nlats)
        eralat = eralats[j]
        colrange = latmap[lat2col(eralat+erares/2, res):lat2col(eralat-erares/2, res)-1]
        for i = 1:nlons
            meanGTI[i,j] == 0 && meanDNI[i,j] == 0 && continue
            GTI = solarGTI[:, i, j]
            DNI = solarDNI[:, i, j]
            eralon = eralons[i]
            # get all high resolution row and column indexes within this ERA5 cell
            rowrange = lonmap[lon2row(eralon-erares/2, res):lon2row(eralon+erares/2, res)-1]

            for c in colrange, r in rowrange
                (c == 0 || r == 0) && continue
                reg = regions[r,c]
                area = cellarea[c]

                class = pvclass[i,j]
                # can't use elseif here, probably some overlap in the masks
                # @views is needed to make sure increment_windCF!() works with matrix slices
                # also faster since it avoids making copies
                @views if reg > 0 && class > 0
                    if mask_rooftop[r,c] > 0
                        # pvrooftopmap[r,c] = meanGTI[i,j]
                        pvrooftopmap[r,c] = class
                    elseif mask_plantA[r,c] > 0
                        # pvmap[r,c] = meanGTI[i,j]
                        pvmap[r,c] = class
                    elseif mask_plantB[r,c] > 0
                        # pvmap[r,c] = class
                    end
                end

                class = cspclass[i,j]
                # @views is needed to make sure increment_windCF!() works with matrix slices
                # also faster since it avoids making copies
                @views if reg > 0 && class > 0
                    if mask_plantA[r,c] > 0
                        # cspmap[r,c] = meanGTI[i,j]
                        cspmap[r,c] = class
                    elseif mask_plantB[r,c] > 0
                        # cspmap[r,c] = class
                    end
                end
            end
        end
        next!(updateprogress)
    end

    return pvmap, pvrooftopmap, cspmap
end

max_stars_repo_name: henrikhodel/GlobalEnergyGIS

max_issues_repo_name: henrikhodel/GlobalEnergyGIS

max_forks_repo_name: henrikhodel/GlobalEnergyGIS


text:
using Serialization

export domain, UnitDomain, project, alignfirst, ArrayDomain

import Base: isempty, getindex, intersect, ==, size, length, ndims

"""
    domain(x::T)

Returns metadata about `x`. This metadata will be in the `domain`
field of a Chunk object when an object of type `T` is created as
the result of evaluating a Thunk.
"""
function domain end

"""
    UnitDomain

Default domain -- has no information about the value
"""
struct UnitDomain end

"""
If no `domain` method is defined on an object, then
we use the `UnitDomain` on it. A `UnitDomain` is indivisible.
"""
domain(x::Any) = UnitDomain()

###### Chunk ######

"""
    Chunk

A reference to a piece of data located on a remote worker. `Chunk`s are
typically created with `Dagger.tochunk(data)`, and the data can then be
accessed from any worker with `collect(::Chunk)`. `Chunk`s are
serialization-safe, and use distributed refcounting (provided by
`MemPool.DRef`) to ensure that the data referenced by a `Chunk` won't be GC'd,
as long as a reference exists on some worker.

Each `Chunk` is associated with a given `Dagger.Processor`, which is (in a
sense) the processor that "owns" or contains the data. Calling
`collect(::Chunk)` will perform data movement and conversions defined by that
processor to safely serialize the data to the calling worker.

## Constructors
See [`tochunk`](@ref).
"""
mutable struct Chunk{T, H, P<:Processor, S<:AbstractScope}
    chunktype::Type{T}
    domain
    handle::H
    processor::P
    scope::S
    persist::Bool
end

domain(c::Chunk) = c.domain
chunktype(c::Chunk) = c.chunktype
persist!(t::Chunk) = (t.persist=true; t)
shouldpersist(p::Chunk) = t.persist
processor(c::Chunk) = c.processor
affinity(c::Chunk) = affinity(c.handle)

function unrelease(c::Chunk{<:Any,DRef})
    # set spilltodisk = true if data is still around
    try
        destroyonevict(c.handle, false)
        return c
    catch err
        if isa(err, KeyError)
            return nothing
        else
            rethrow(err)
        end
    end
end
unrelease(c::Chunk) = c

Base.:(==)(c1::Chunk, c2::Chunk) = c1.handle == c2.handle
Base.hash(c::Chunk, x::UInt64) = hash(c.handle, x)

collect_remote(chunk::Chunk) =
    move(chunk.processor, OSProc(), poolget(chunk.handle))

function collect(ctx::Context, chunk::Chunk; options=nothing)
    # delegate fetching to handle by default.
    if chunk.handle isa DRef && !(chunk.processor isa OSProc)
        return remotecall_fetch(collect_remote, chunk.handle.owner, chunk)
    elseif chunk.handle isa FileRef
        return poolget(chunk.handle)
    else
        return move(chunk.processor, OSProc(), chunk.handle)
    end
end
collect(ctx::Context, ref::DRef; options=nothing) =
    move(OSProc(ref.owner), OSProc(), ref)
collect(ctx::Context, ref::FileRef; options=nothing) =
    poolget(ref) # FIXME: Do move call
function Base.fetch(chunk::Chunk; raw=false)
    if raw
        poolget(chunk.handle)
    else
        collect(chunk)
    end
end

# Unwrap Chunk, DRef, and FileRef by default
move(from_proc::Processor, to_proc::Processor, x::Chunk) =
    move(from_proc, to_proc, x.handle)
move(from_proc::Processor, to_proc::Processor, x::Union{DRef,FileRef}) =
    move(from_proc, to_proc, poolget(x))

# Determine from_proc when unspecified
move(to_proc::Processor, chunk::Chunk) =
    move(chunk.processor, to_proc, chunk)
move(to_proc::Processor, d::DRef) =
    move(OSProc(d.owner), to_proc, d)
move(to_proc::Processor, x) =
    move(OSProc(), to_proc, x)

### ChunkIO
affinity(r::DRef) = OSProc(r.owner)=>r.size
# this previously returned a vector with all machines that had the file cached
# but now only returns the owner and size, for consistency with affinity(::DRef),
# see #295
affinity(r::FileRef) = OSProc(1)=>r.size

### Mutation

"Wraps `x` in a `Chunk` on `proc`, scoped to `scope`, which allows `x` to be mutated by tasks that use it."
macro mutable(proc, scope, x)
    :(Dagger.tochunk($(esc(x)), $(esc(proc)), $(esc(scope))))
end
"Creates a mutable `Chunk` on `proc`, scoped to exactly `proc`."
macro mutable(proc, x)
    quote
        let proc = $(esc(proc))
            let scope = proc isa OSProc ? Dagger.ProcessScope(proc.pid) : Dagger.ExactScope(proc)
                Dagger.@mutable proc scope $(esc(x))
            end
        end
    end
end
"Creates a mutable `Chunk` on the current worker."
macro mutable(x)
    :(Dagger.@mutable OSProc() Dagger.ProcessScope() $(esc(x)))
end

"""
Maps a value to one of multiple distributed "mirror" values automatically when
used as a thunk argument. Construct using `@shard` or `shard`.
"""
struct Shard
    chunks::Dict{Processor,Chunk}
end

"""
    shard(f; kwargs...) -> Chunk{Shard}

Executes `f` on all workers in `workers`, wrapping the result in a
process-scoped `Chunk`, and constructs a `Chunk{Shard}` containing all of these
`Chunk`s on the current worker.

Keyword arguments:
- `procs` -- The list of processors to create pieces on. May be any iterable container of `Processor`s.
- `workers` -- The list of workers to create pieces on. May be any iterable container of `Integer`s.
- `per_thread::Bool=false` -- If `true`, creates a piece per each thread, rather than a piece per each worker.
"""
function shard(f; procs=nothing, workers=nothing, per_thread=false)
    if procs === nothing
        if workers !== nothing
            procs = [OSProc(w) for w in workers]
        else
            procs = lock(Sch.eager_context()) do
                copy(Sch.eager_context().procs)
            end
        end
        if per_thread
            _procs = ThreadProc[]
            for p in procs
                append!(_procs, filter(p->p isa ThreadProc, get_processors(p)))
            end
            procs = _procs
        end
    else
        if workers !== nothing
            throw(ArgumentError("Cannot combine `procs` and `workers`"))
        elseif per_thread
            throw(ArgumentError("Cannot combine `procs` and `per_thread=true`"))
        end
    end
    isempty(procs) && throw(ArgumentError("Cannot create empty Shard"))
    scopes = [p isa OSProc ? ProcessScope(p) : ExactScope(p) for p in procs]
    thunks = [proc=>Dagger.@spawn single=Dagger.get_parent(proc).pid _shard_inner(f, proc, scope) for (proc,scope) in zip(procs,scopes)]
    shard = Shard(Dict{Processor,Chunk}(thunk[1]=>fetch(thunk[2])[] for thunk in thunks))
    scope = UnionScope(scopes)
    Dagger.tochunk(shard, OSProc(), scope)
end
function _shard_inner(f, proc, scope)
    Ref(Dagger.@mutable proc scope f())
end

"Creates a `Shard`. See [`Dagger.shard`](@ref) for details."
macro shard(exs...)
    opts = esc.(exs[1:end-1])
    ex = exs[end]
    quote
        let f = ()->$(esc(ex))
            $shard(f; $(opts...))
        end
    end
end

function move(from_proc::Processor, to_proc::Processor, shard::Shard)
    # Match either this proc or some ancestor
    # N.B. This behavior may bypass the piece's scope restriction
    proc = to_proc
    if haskey(shard.chunks, proc)
        return shard.chunks[proc]
    end
    parent = Dagger.get_parent(proc)
    while parent != proc
        proc = parent
        parent = Dagger.get_parent(proc)
        if haskey(shard.chunks, proc)
            return shard.chunks[proc]
        end
    end

    throw(KeyError(to_proc))
end
function move(from_proc::Processor, to_proc::Processor, x::Chunk{Shard})
    piece = remotecall_fetch(x.handle.owner, x.handle, from_proc, to_proc) do ref, from_proc, to_proc
        shard = MemPool.poolget(ref)
        move(from_proc, to_proc, shard)
    end::Chunk
    move(from_proc, to_proc, piece)
end
Base.map(f, cs::Chunk{Shard}) = map(f, fetch(cs; raw=true))
Base.map(f, s::Shard) = [Dagger.spawn(f, c) for c in values(s.chunks)]

### Core Stuff

"""
    tochunk(x, proc; persist=false, cache=false) -> Chunk

Create a chunk from sequential object `x` which resides on `proc`.
"""
function tochunk(x::X, proc::P=OSProc(), scope::S=AnyScope(); persist=false, cache=false) where {X,P,S}
    ref = poolset(x, destroyonevict=persist ? false : cache)
    Chunk{X,typeof(ref),P,S}(X, domain(x), ref, proc, scope, persist)
end
tochunk(x::Union{Chunk, Thunk}, proc=nothing, scope=nothing) = x

# Check to see if the node is set to persist
# if it is foce can override it
function free!(s::Chunk{X,DRef,P,S}; force=true, cache=false) where {X,P,S}
    if force || !s.persist
        if cache
            try
                destroyonevict(s.handle, true) # keep around, but remove when evicted
            catch err
                isa(err, KeyError) || rethrow(err)
            end
        end
    end
end
free!(x; force=true,cache=false) = x # catch-all for non-chunks

function savechunk(data, dir, f)
    sz = open(joinpath(dir, f), "w") do io
        serialize(io, MemPool.MMWrap(data))
        return position(io)
    end
    fr = FileRef(f, sz)
    proc = OSProc()
    scope = AnyScope() # FIXME: Scoped to this node
    Chunk{typeof(data),typeof(fr),typeof(proc),typeof(scope)}(typeof(data), domain(data), fr, proc, scope, true)
end

Base.@deprecate_binding AbstractPart Union{Chunk, Thunk}
Base.@deprecate_binding Part Chunk
Base.@deprecate parts(args...) chunks(args...)
Base.@deprecate part(args...) tochunk(args...)
Base.@deprecate parttype(args...) chunktype(args...)


text:
"""
  LsqProblem(residuals)

An LsqProblem ([Non-linear] Least squares Problem), is used to represent the
mathematical problem of finding the smallest possible sum of squared values of
the elements of the residual function. The problem is defined by `residual` which is
an appropriate objective type (for example `NonDiffed`, `OnceDiffed`, ...) for the
types of algorithm to be used.

Options are stored in `options` and are of the `NLsqOptions` type. See more information
about options using `?NLsqOptions`.

The package NLSolversAD.jl adds automatic conversion of problems to match algorithms
that require higher order derivates than provided by the user. It also adds AD
constructors for a target number of derivatives.
"""
struct LsqProblem{R<:ObjWrapper, B, M, C}
  residuals::R
  bounds::B
  manifold::M
  constraints::C
end

struct LsqOptions{Tmi}
  maxiter::Tmi
end


text:
using Requires
@init @require CUDAnative = "be33ccc6-a3ff-5ff2-a52e-74243cff1e17" begin
    using .CUDAnative
end

@kernel function linearcombination!(Q, cs, Xs, increment::Bool)
    i = @index(Global, Linear)
    if !increment
        @inbounds Q[i] = -zero(eltype(Q))
    end
    @inbounds for j in 1:length(cs)
        Q[i] += cs[j] * Xs[j][i]
    end
end

max_stars_repo_name: yairchn/CLIMA

max_issues_repo_name: yairchn/CLIMA

max_forks_repo_name: yairchn/CLIMA


text:
export process, process_ind, process_f

function process_ind(ind::CGPInd, inputs::AbstractArray)
    # Set inputs
    @inbounds for i in eachindex(inputs)
        ind.buffer[i] = inputs[i]
    end
    # Process
    @inbounds for i in eachindex(ind.nodes)
        if ind.nodes[i].active
            ind.buffer[i] = ind.nodes[i].f(
                ind.buffer[ind.nodes[i].x],
                ind.buffer[ind.nodes[i].y],
                ind.nodes[i].p
            )
        end
    end
    # Get outputs
    [ind.buffer[i] for i in ind.outputs]
end

"""
    function process(
        encoder::CGPInd,
        reducer::AbstractReducer,
        controller::CGPInd,
        controller_config::NamedTuple,
        inp::AbstractArray
    )

Process function chaining encoder, features projection and controller.
Both the encoder and the controller are CGP individuals.
"""
function process(
    encoder::CGPInd,
    reducer::AbstractReducer,
    controller::CGPInd,
    controller_config::NamedTuple,
    inp::AbstractArray
)
    out = CartesianGeneticProgramming.process(encoder, inp)
    features = reducer.reduct(out, reducer.parameters)
    features_flatten = collect(Iterators.flatten(Iterators.flatten(features)))
    if controller_config.n_cst_inputs > 1
        push!(features_flatten, 0:1.0/(controller_config.n_cst_inputs-1):1.0...)
    end
    CartesianGeneticProgramming.process(controller, features_flatten)
end

"""
    function process_f(
        encoder::CGPInd,
        reducer::AbstractReducer,
        controller::CGPInd,
        controller_config::NamedTuple,
        inp::AbstractArray
    )

Process function chaining encoder, features projection and controller.
Both the encoder and the controller are CGP individuals.
Return both the created feature vector and the output.
"""
function process_f(
    encoder::CGPInd,
    reducer::AbstractReducer,
    controller::CGPInd,
    controller_config::NamedTuple,
    inp::AbstractArray
)
    out = CartesianGeneticProgramming.process(encoder, inp)
    features = reducer.reduct(out, reducer.parameters)
    # features_flatten = collect(Iterators.flatten(features))
    features_flatten = collect(Iterators.flatten(Iterators.flatten(features)))
    if controller_config.n_cst_inputs > 1
        push!(features_flatten, 0:1.0/(controller_config.n_cst_inputs-1):1.0...)
    end
    return features, CartesianGeneticProgramming.process(controller, features_flatten)
end

"""
    function process_full(
        encoder::CGPInd,
        reducer::AbstractReducer,
        controller::CGPInd,
        controller_config::NamedTuple,
        inp::AbstractArray
    )

Process function chaining encoder, features projection and controller.
Both the encoder and the controller are CGP individuals.
Return the encoder output, the created feature vector and the output.
"""
function process_full(
    encoder::CGPInd,
    reducer::AbstractReducer,
    controller::CGPInd,
    controller_config::NamedTuple,
    inp::AbstractArray
)
    enco_out = CartesianGeneticProgramming.process(encoder, inp)
    features = reducer.reduct(enco_out, reducer.parameters)
    # features_flatten = collect(Iterators.flatten(features))
    features_flatten = collect(Iterators.flatten(Iterators.flatten(features)))
    if controller_config.n_cst_inputs > 1
        push!(features_flatten, 0:1.0/(controller_config.n_cst_inputs-1):1.0...)
    end
    cont_out = CartesianGeneticProgramming.process(controller, features_flatten)
    return enco_out, features, cont_out
end

function process(
    reducer::AbstractReducer,
    controller::CGPInd,
    inp::AbstractArray,
)
    features = reducer.reduct(inp, reducer.parameters)
    features_flatten = collect(Iterators.flatten(Iterators.flatten(features)))
    CartesianGeneticProgramming.process(controller, features_flatten)
end


text:
using JSON, Sockets, Pkg, SymbolServer

init_request = """
{
    "jsonrpc":"2.0",
    "id":0,
    "method":"initialize",
    "params":{"processId":9902,
              "rootPath":null,
              "rootUri":null,
              "capabilities":{"workspace":{"applyEdit":true,"workspaceEdit":{"documentChanges":true},"didChangeConfiguration":{"dynamicRegistration":false},"didChangeWatchedFiles":{"dynamicRegistration":false},"symbol":{"dynamicRegistration":true},"executeCommand":{"dynamicRegistration":true}},"textDocument":{"synchronization":{"dynamicRegistration":true,"willSave":true,"willSaveWaitUntil":true,"didSave":true},"completion":{"dynamicRegistration":true,"completionItem":{"snippetSupport":true}},"hover":{"dynamicRegistration":true},"signatureHelp":{"dynamicRegistration":true},"references":{"dynamicRegistration":true},"documentHighlight":{"dynamicRegistration":true},"documentSymbol":{"dynamicRegistration":true},"formatting":{"dynamicRegistration":true},"rangeFormatting":{"dynamicRegistration":true},"onTypeFormatting":{"dynamicRegistration":true},"definition":{"dynamicRegistration":true},"codeAction":{"dynamicRegistration":true},"codeLens":{"dynamicRegistration":true},"documentLink":{"dynamicRegistration":true},"rename":{"dynamicRegistration":true}}},
              "trace":"off"}
}
"""

init_response_json = JSON.parse("""
{
    "id": 0,
    "jsonrpc": "2.0",
    "result": {
        "capabilities": {
            "textDocumentSync": 2,
            "hoverProvider": true,
            "completionProvider": {
                "resolveProvider": false,
                "triggerCharacters": [
                    "."
                ]
            },
            "signatureHelpProvider": {
                "triggerCharacters": [
                    "("
                ]
            },
            "definitionProvider": true,
            "typeDefinitionProvider": false,
            "implementationProvider": false,
            "referencesProvider": true,
            "documentHighlightProvider": false,
            "documentSymbolProvider": true,
            "workspaceSymbolProvider": true,
            "codeActionProvider": true,
            "documentFormattingProvider": true,
            "documentRangeFormattingProvider": false,
            "renameProvider": true,
            "documentLinkProvider": {
                "resolveProvider": false
            },
            "colorProvider": false,
            "executeCommandProvider": {
                "commands": []
            },
            "workspace": {
                "workspaceFolders": {
                    "supported": true,
                    "changeNotifications": true
                }
            },
            "experimental": null
        }
    }
}
""")

if Sys.iswindows()
    global_socket_name = "\\\\.\\pipe\\julia-language-server-testrun"
elseif Sys.isunix()
    global_socket_name = joinpath(tempdir(), "julia-language-server-testrun")
else
    error("Unknown operating system.")
end

@async begin
    server = listen(global_socket_name)
    try
        sock = accept(server)
        try
            ls = LanguageServerInstance(sock, sock, false, dirname(Pkg.Types.Context().env.project_file), first(Base.DEPOT_PATH), Dict())
            try
                run(ls)
            catch err
                Base.display_error(stderr, err, catch_backtrace())
                rethrow()
            end
        finally
            close(sock)
        end
    finally
        close(server)
    end
end

sleep(1)

client = connect(global_socket_name)
try
    LanguageServer.write_transport_layer(client, init_request)

    msg = LanguageServer.read_transport_layer(client)

    msg_json = JSON.parse(msg)

    @test_broken init_response_json == msg_json
    @test msg_json["result"]["capabilities"]["typeDefinitionProvider"] == false
    @test msg_json["result"]["capabilities"]["renameProvider"] == true
finally
    close(client)
end


text:
# All types `J <: AbstractJob` must implement:
#
# - `J(submission::JobSubmission)`: create an instance of `J` from `submission`
# - `submission(job::J)`: return the `JobSubmission` used to create `job`
# - `isvalid(submission::JobSubmission, ::Type{J})`: return true if `submission` is a valid input for `J`, false otherwise
# - `Base.run(job::J)`: execute `job`
# - `Base.summary(job::J)`: a short string identifying/describing `job`

abstract type AbstractJob end

reply_status(job::AbstractJob, args...; kwargs...) = reply_status(submission(job), "Nanosoldier/$(nameof(typeof(job)))", args...; kwargs...)
reply_comment(job::AbstractJob, args...; kwargs...) = reply_comment(submission(job), args...; kwargs...)
upload_report_repo!(job::AbstractJob, args...; kwargs...) = upload_report_repo!(submission(job), args...; kwargs...)

function branchref(config::Config, reponame::AbstractString, branchname::AbstractString)
    shastr = GitHub.branch(reponame, branchname; auth=config.auth).commit.sha
    return BuildRef(reponame, shastr)
end

function tagref(config::Config, reponame::AbstractString, tagname::AbstractString)
    shastr = GitHub.tag(reponame, tagname; auth=config.auth).object["sha"]
    return BuildRef(reponame, shastr)
end

datedirname(date::Dates.Date) = joinpath(Dates.format(date, dateformat"yyyy-mm"),
                                         Dates.format(date, dateformat"dd"))

include("BenchmarkJob.jl")
include("PkgEvalJob.jl")


text:
module SimLib

using Dates
import JLD2
using LinearAlgebra
using Reexport
using Statistics
using XXZNumerics

export logmsg, path_prefix, parse_geometry, geometry_from_density, meandrop, stddrop
export SaveLocation, descriptor, datapath, create, save, load, load_or_create

# simplify type definitions
const FArray{N} = Array{Float64, N} where N
const Maybe{T} = Union{Missing, T} where T

include("general.jl")
include("data.jl")
include("positions.jl")
include("exact_diagonalization/exact_diagonalization.jl")
include("eigenstate_occupation.jl")
include("eigenstate_locality.jl")
include("levels.jl")
include("operator_diagonal.jl")
include("ensembles.jl")
include("levelspacingratio.jl")
include("halfchainentropy_zblock.jl")
include("inverseparticipationratio.jl")

@reexport using .Positions
@reexport using .ED
@reexport using .Levels
@reexport using .EON
@reexport using .OPDiag
@reexport using .LSR
@reexport using .EL
@reexport using .Ensembles
@reexport using .HCE_ZBlock_Module
@reexport using .IPR

end


text:
# MIT License

# Copyright (c) Microsoft Corporation.

# Permission is hereby granted, free of charge, to any person obtaining a copy
# of this software and associated documentation files (the "Software"), to deal
# in the Software without restriction, including without limitation the rights
# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
# copies of the Software, and to permit persons to whom the Software is
# furnished to do so, subject to the following conditions:

# The above copyright notice and this permission notice shall be included in all
# copies or substantial portions of the Software.

# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
# SOFTWARE

"""Abstract type representing any evaluated CSG structure."""
abstract type CSGTree{T} <: Primitive{T} end
export CSGTree

"""
    ComplementNode{T,C<:CSGTree{T}} <: CSGTree{T}

An evaluated complement node within the CSG tree, must be the second child of a [`IntersectionNode`](@ref) forming a subtraction.
"""
struct ComplementNode{T,C<:CSGTree{T}} <: CSGTree{T}
    child::C

    function ComplementNode(child::C) where {T<:Real,C<:CSGTree{T}}
        return new{T,C}(child)
    end
end
Base.show(io::IO, a::ComplementNode{T}) where {T} = print(io, "Complement($(a.child))")
BoundingBox(a::ComplementNode{T}) where {T<:Real} = BoundingBox(a.child)

"""
    UnionNode{T,L<:CSGTree{T},R<:CSGTree{T}} <: CSGTree{T}

An evaluated union node within the CSG tree.
"""
struct UnionNode{T,L<:CSGTree{T},R<:CSGTree{T}} <: CSGTree{T}
    leftchild::L
    rightchild::R
    bbox::BoundingBox{T}

    function UnionNode(a::L, b::R) where {T<:Real,L<:CSGTree{T},R<:CSGTree{T}}
        # union should never contain a complement so should be fine
        return new{T,L,R}(a, b, union(BoundingBox(a), BoundingBox(b)))
    end
end
Base.show(io::IO, a::UnionNode{T}) where {T} = print(io, "Union($(a.leftchild), $(a.rightchild))")

"""
    IntersectionNode{T,L<:CSGTree{T},R<:CSGTree{T}} <: CSGTree{T}

An evaluated intersection node within the CSG tree.
"""
struct IntersectionNode{T,L<:CSGTree{T},R<:CSGTree{T}} <: CSGTree{T}
    leftchild::L
    rightchild::R
    bbox::BoundingBox{T}

    function IntersectionNode(a::L, b::R) where {T<:Real,L<:CSGTree{T},R<:CSGTree{T}}
        # normal intersection is fine for most nodes
        return new{T,L,R}(a, b, intersection(BoundingBox(a), BoundingBox(b)))
    end

    function IntersectionNode(a::L, b::R) where {T<:Real,L<:CSGTree{T},R<:ComplementNode{T}}
        # this is a subtraction so just take the original unclipped bounding box for simplicity
        return new{T,L,R}(a, b, BoundingBox(a))
    end
end
Base.show(io::IO, a::IntersectionNode{T}) where {T} = print(io, "Intersection($(a.leftchild), $(a.rightchild))")

"""
    LeafNode{T,S<:ParametricSurface{T}} <: CSGTree{T}

An evaluated leaf node in the CSG tree, `geometry` attribute which contains a [`ParametricSurface`](@ref) of type `S`.
The leaf node also has a transform associated which is the composition of all nodes above it in the tree.
As such, transforming points from the geometry using this transform puts them in world space, and transforming rays by the inverse transform puts them in object space.
"""
struct LeafNode{T,S<:ParametricSurface{T,3}} <: CSGTree{T}
    geometry::S
    transform::RigidBodyTransform{T}
    invtransform::RigidBodyTransform{T}
    bbox::BoundingBox{T}

    function LeafNode(a::S, transform::RigidBodyTransform{T}) where {T<:Real,S<:ParametricSurface{T,3}}
        # store the transformed bounding box so nodes higher in the tree have correct global space bounding boxes
        return new{T,S}(a, transform, inv(transform), BoundingBox(a, transform))
    end
end
Base.show(io::IO, a::LeafNode{T}) where {T} = a.transform != identitytransform(T) ? print(io, "Leaf($(a.geometry), $(a.transform))") : print(io, "Leaf($(a.geometry))")

BoundingBox(a::CSGTree{T}) where {T<:Real} = a.bbox

"""
    CSGGenerator{T<:Real}

This is the type you should use when making CSG objects.
This type allows for the construction of [`CSGTree`](@ref) objects with different transforms.
When the generator is evaluated, all transforms are propagated down to the [`LeafNode`](@ref)s and stored there.

# Example
```julia
a = Cylinder(1.0,1.0)
b = Plane([0.0,0.0,1.0], [0.0,0.0,0.0])
generator = csgintersection(a,b)
# now make a csg object that can be ray traced
csgobj = generator(RigidBodyTransform(1.0,1.0,2.0))
```
"""
struct CSGGenerator{T<:Real}
    f::Function
end
export CSGGenerator

function (a::CSGGenerator{T})(transform::RigidBodyTransform{T})::CSGTree{T} where {T<:Real}
    a.f(transform)
end
function (a::CSGGenerator{T})()::CSGTree{T} where {T<:Real}
    a.f(identitytransform(T))
end

"""
    leaf(surf::ParametricSurface{T}, transform::RigidBodyTransform{T} = identitytransform(T)) -> CSGGenerator{T}

Create a leaf node from a parametric surface with a given transform.
"""
function leaf(surf::S, transform::RigidBodyTransform{T} = identitytransform(T))::CSGGenerator{T} where {T<:Real,S<:ParametricSurface{T}}
    return CSGGenerator{T}((parenttransform) -> LeafNode(surf, parenttransform * transform))
end
"""
    leaf(surf::CSGGenerator{T}, transform::RigidBodyTransform{T} = identitytransform(T)) -> CSGGenerator{T}

Create a (pseudo) leaf node from another CSGGenerator, this is useful if you want multiple copies of a premade CSG structure with different transforms, for example in an MLA.
"""
function leaf(n::CSGGenerator{T}, transform::RigidBodyTransform{T} = identitytransform(T))::CSGGenerator{T} where {T<:Real}
    return CSGGenerator{T}((parenttransform) -> n(parenttransform * transform))
end
export leaf

# CSG objects are trees, not graphs, since it makes no sense to reuse intermediate results. Hence no need to keep track of common subexpressions.
# Static arrays are much faster than mutable arrays so want node transforms to be static. Unfortunately the transformations cascade from the root to the leaves but the nodes have to be created from the leaves to the roots.
# Wrap the csg operations in function that delays evaluation until the transform has been computed.

"""
    csgintersection(a::CSGGenerator{T} b::CSGGenerator{T}, transform::RigidBodyTransform{T} = identitytransform(T)) -> CSGGenerator{T}

Create a binary node in the CSG tree representing an intersection between a and b.
A shortcut method for `a` and `b` as [`ParametricSurface`](@ref)s is also available.

![Intersect Image](https://upload.wikimedia.org/wikipedia/commons/0/0b/Boolean_intersect.PNG)
"""
csgintersection(a::CSGGenerator{T}, b::CSGGenerator{T}, transform::RigidBodyTransform{T} = identitytransform(T)) where {T<:Real} = CSGGenerator{T}((parenttransform) -> IntersectionNode(a(parenttransform * transform), b(parenttransform * transform)))
csgintersection(a::ParametricSurface{T}, b::ParametricSurface{T}, transform::RigidBodyTransform{T} = identitytransform(T)) where {T<:Real} = csgintersection(leaf(a), leaf(b), transform)
csgintersection(a::ParametricSurface{T}, b::CSGGenerator{T}, transform::RigidBodyTransform{T} = identitytransform(T)) where {T<:Real} = csgintersection(leaf(a), b, transform)
csgintersection(a::CSGGenerator{T}, b::ParametricSurface{T}, transform::RigidBodyTransform{T} = identitytransform(T)) where {T<:Real} = csgintersection(a, leaf(b), transform)

export csgintersection

"""
    csgunion(a::CSGGenerator{T}, b::CSGGenerator{T}, transform::RigidBodyTransform{T} = identitytransform(T)) -> CSGGenerator{T}

Create a binary node in the CSG tree representing a union between a and b.
A shortcut method for `a` and `b` as [`ParametricSurface`](@ref)s is also available.

![Union Image](https://upload.wikimedia.org/wikipedia/commons/4/4a/Boolean_union.PNG)
"""
csgunion(a::CSGGenerator{T}, b::CSGGenerator{T}, transform::RigidBodyTransform{T} = identitytransform(T)) where {T<:Real} = CSGGenerator{T}((parenttransform) -> UnionNode(a(parenttransform * transform), b(parenttransform * transform)))
csgunion(a::ParametricSurface{T}, b::ParametricSurface{T}, transform::RigidBodyTransform{T} = identitytransform(T)) where {T<:Real} = csgunion(leaf(a), leaf(b), transform)
csgunion(a::ParametricSurface{T}, b::CSGGenerator{T}, transform::RigidBodyTransform{T} = identitytransform(T)) where {T<:Real} = ccsgunion(leaf(a), b, transform)
csgunion(a::CSGGenerator{T}, b::ParametricSurface{T}, transform::RigidBodyTransform{T} = identitytransform(T)) where {T<:Real} = csgunion(a, leaf(b), transform)

export csgunion

"""
    csgdifference(a::CSGGenerator{T}, b::CSGGenerator{T}, transform::RigidBodyTransform{T} = identitytransform(T)) -> CSGGenerator{T}

Create a binary node in the CSG tree representing the difference of a and b, essentially a - b.
A shortcut method for `a` and `b` as [`ParametricSurface`](@ref)s is also available.

![Difference Image](https://upload.wikimedia.org/wikipedia/commons/8/86/Boolean_difference.PNG)
"""
csgdifference(a::CSGGenerator{T}, b::CSGGenerator{T}, transform::RigidBodyTransform{T} = identitytransform(T)) where {T<:Real} = CSGGenerator{T}((parenttransform) -> IntersectionNode(a(parenttransform * transform), ComplementNode(b(parenttransform * transform))))
csgdifference(a::CSGGenerator{T}, b::ParametricSurface{T}, transform::RigidBodyTransform{T} = identitytransform(T)) where {T<:Real} = csgdifference(a, leaf(b), transform)
csgdifference(a::ParametricSurface{T}, b::CSGGenerator{T}, transform::RigidBodyTransform{T} = identitytransform(T)) where {T<:Real} = csgdifference(leaf(a), b, transform)
csgdifference(a::ParametricSurface{T}, b::ParametricSurface{T}, transform::RigidBodyTransform{T} = identitytransform(T)) where {T<:Real} = csgdifference(leaf(a), leaf(b), transform)

export csgdifference

function evalcsg(a::UnionNode{T}, ray::AbstractRay{T,N}, normalreverse::Bool = false)::Union{EmptyInterval{T},DisjointUnion{T},Interval{T}} where {T<:Real,N}
    if !doesintersect(a.bbox, ray)
        return EmptyInterval(T)
    end
    l = evalcsg(a.leftchild, ray, normalreverse)
    if l isa Interval{T}
        if lower(l) isa RayOrigin{T} && upper(l) isa Infinity{T}
            return l
        end
    end
    r = evalcsg(a.rightchild, ray, normalreverse)
    if l isa EmptyInterval{T}
        return r
    end
    if r isa Interval{T}
        if lower(r) isa RayOrigin{T} && upper(r) isa Infinity{T}
            return r
        end
    end
    if r isa EmptyInterval{T}
        return l
    end
    return intervalunion(l, r)
end

function evalcsg(a::IntersectionNode{T}, ray::AbstractRay{T,N}, normalreverse::Bool = false)::Union{EmptyInterval{T},DisjointUnion{T},Interval{T}} where {T<:Real,N}
    if !doesintersect(a.bbox, ray)
        return EmptyInterval(T)
    end
    l = evalcsg(a.leftchild, ray, normalreverse)
    if l isa EmptyInterval{T}
        return EmptyInterval(T)
    else
        r = evalcsg(a.rightchild, ray, normalreverse)
        if r isa EmptyInterval{T}
            return EmptyInterval(T)
        else
            return intervalintersection(l, r)
        end
    end
end

function evalcsg(a::ComplementNode{T}, ray::AbstractRay{T,N}, normalreverse::Bool = false)::Union{EmptyInterval{T},DisjointUnion{T},Interval{T}} where {T<:Real,N}
    return intervalcomplement(evalcsg(a.child, ray, !normalreverse))
end

function evalcsg(a::LeafNode{T}, ray::AbstractRay{T,N}, normalreverse::Bool = false)::Union{EmptyInterval{T},DisjointUnion{T},Interval{T}} where {T<:Real,N}
    # the bounding box is in global coordinates so use the un-transformed ray
    if !doesintersect(a.bbox, ray)
        return EmptyInterval(T)
    end
    intscts = surfaceintersection(a.geometry, a.invtransform * ray)
    if !(intscts isa EmptyInterval{T})
        if normalreverse
            intscts = reversenormal(intscts)
        end
        if intscts isa Interval{T}
            return a.transform * intscts
        elseif intscts isa DisjointUnion{T}
            for i in eachindex(intervals(intscts))
                intervals(intscts)[i] = a.transform * intervals(intscts)[i]
            end
            return intscts
        else
            throw(ErrorException("EvalCSG error - returned type not an interval or disjoint union"))
        end
    else
        return EmptyInterval(T)
    end
end

"""
    surfaceintersection(obj::CSGTree{T}, r::AbstractRay{T,N})

Calculates the intersection of `r` with CSG object, `obj`.

Returns an [`EmptyInterval`](@ref) if there is no intersection, an [`Interval`](@ref) if there is one or two interesections and a [`DisjointUnion`](@ref) if there are more than two intersections.

The ray is intersected with the [`LeafNode`](@ref)s that make up the CSG object and the resulting `Interval`s and `DisjointUnion`s are composed with the same boolean operations to give a final result.
The ray is transformed by the inverse of the transform associated with the leaf node to put it in _object space_ for that node before the intersection is carried out, typically this _object space_ is centered at the origin, but may differ for each primitive.

Some intersections are culled without actually evaluating them by first checking if the ray intersects the [`BoundingBox`](@ref) of each node in the [`CSGTree`](@ref), this can substantially improve performance in some cases.
"""
surfaceintersection(element::CSGTree{T}, r::AbstractRay{T,N}) where {T<:Real,N} = evalcsg(element, r)

point(::EmptyInterval) = nothing
normal(::EmptyInterval) = nothing
point(a::DisjointUnion) = point(closestintersection(a))
normal(a::DisjointUnion) = normal(closestintersection(a))
point(a::Interval) = point(closestintersection(a))
normal(a::Interval) = point(closestintersection(a))

"""
    onsurface(obj::CSGTree{T}, point::SVector{3,T}) -> Bool
    onsurface(obj::CSGTree{T}, x::T, y::T, z::T) -> Bool

Tests whether a 3D point in world space is _on_ the surface (i.e. shell) of `obj`.
"""
onsurface(a::CSGTree{T}, x::T, y::T, z::T) where {T<:Real} = onsurface(a, SVector{3,T}(x, y, z))
onsurface(a::ComplementNode{T}, point::SVector{3,T}) where {T<:Real} = onsurface(a.child, point)
onsurface(a::IntersectionNode{T}, point::SVector{3,T}) where {T<:Real} = (onsurface(a.leftchild, point) && inside(a.rightchild, point)) || (onsurface(a.rightchild, point) && inside(a.leftchild, point))
onsurface(a::UnionNode{T}, point::SVector{3,T}) where {T<:Real} = (onsurface(a.leftchild, point) && !inside(a.rightchild, point)) || (onsurface(a.rightchild, point) && !inside(a.leftchild, point))
onsurface(a::LeafNode{T}, point::SVector{3,T}) where {T<:Real} = onsurface(a.geometry, a.invtransform * point)

"""
    inside(obj::CSGTree{T}, point::SVector{3,T}) -> Bool
    inside(obj::CSGTree{T}, x::T, y::T, z::T) -> Bool

Tests whether a 3D point in world space is _inside_ `obj`.
"""
inside(a::CSGTree{T}, x::T, y::T, z::T) where {T<:Real} = inside(a, SVector{3,T}(x, y, z))
inside(a::ComplementNode{T}, point::SVector{3,T}) where {T<:Real} = !inside(a.child, point)
inside(a::IntersectionNode{T}, point::SVector{3,T}) where {T<:Real} = inside(a.bbox, point) && (inside(a.leftchild, point) && inside(a.rightchild, point))
inside(a::UnionNode{T}, point::SVector{3,T}) where {T<:Real} = inside(a.bbox, point) && (inside(a.leftchild, point) || inside(a.rightchild, point))
inside(a::LeafNode{T}, point::SVector{3,T}) where {T<:Real} = inside(a.bbox, point) && inside(a.geometry, a.invtransform * point)

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

struct TrianglePool{T<:Real}
    allocated::Vector{Vector{Triangle{T}}}
    unallocated::Vector{Vector{Triangle{T}}}

    TrianglePool{T}() where {T<:Real} = new{T}(Vector{Vector{Triangle{T}}}(), Vector{Vector{Triangle{T}}}())
end

function allocate!(a::TrianglePool{T}) where {T<:Real}
    if length(a.unallocated) === 0
        push!(a.unallocated, Vector{Triangle{T}}()) #this will extend the array with a new Vector of Triangles.
    end

    temp = pop!(a.unallocated)
    Base.empty!(temp) #resets count in array but doesn't reclaim space
    push!(a.allocated, temp)
    return temp
end

function empty!(a::TrianglePool{T}) where {T<:Real}
    for i in 1:length(a.allocated)
        push!(a.unallocated, pop!(a.allocated))
    end
end

const threadedtrianglepool = [Dict{DataType,TrianglePool}([Float64 => TrianglePool{Float64}()]) for _ in 1:Threads.nthreads()]

function newintrianglepool!(::Type{T} = Float64, tid::Int = Threads.threadid())::Vector{Triangle{T}} where {T<:Real}
    if T ∉ keys(threadedtrianglepool[tid])
        # if the type of the interval pool has changed then we need to refill it with the correct type
        threadedtrianglepool[tid][T] = TrianglePool{T}()
    end
    return allocate!(threadedtrianglepool[tid][T])
end

function emptytrianglepool!(::Type{T} = Float64, tid::Int = Threads.threadid()) where {T}
    if T ∈ keys(threadedtrianglepool[tid])
        OpticSim.empty!(threadedtrianglepool[tid][T])
    end
end

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

# currrently we assume that no triangles span 2 sides of a shape, i.e. any triangle only crosses one shape boundary and so splitting is trivial
function uniontri!(csg::CSGTree{T}, tri::Triangle{T}, triangles::Vector{Triangle{T}}, thisforcoplanar = false) where {T<:Real}
    v1 = vertex(tri, 1)
    v2 = vertex(tri, 2)
    v3 = vertex(tri, 3)
    in1 = inside(csg, v1)
    in2 = inside(csg, v2)
    in3 = inside(csg, v3)
    # second condition for if the vertex is in the plane of the other object, in which case we want to keep one copy of the two coplanar surfaces
    if !thisforcoplanar
        # if we are NOT using this surface for coplanar faces then we want to count the verts on the surface as
        # being inside too in case we have coplanar faces
        in1 |= onsurface(csg, v1)
        in2 |= onsurface(csg, v2)
        in3 |= onsurface(csg, v3)
    end
    if !(in1 || in2 || in3) || (thisforcoplanar && onsurface(csg, v1) && onsurface(csg, v2) && onsurface(csg, v3))
        # whole triangle is valid
        push!(triangles, tri)
    elseif !in1 && !in2 && in3
        splittri2out!(csg, v1, v2, v3, triangles)
    elseif !in2 && !in3 && in1
        splittri2out!(csg, v2, v3, v1, triangles)
    elseif !in1 && !in3 && in2
        splittri2out!(csg, v3, v1, v2, triangles)
    elseif !in1 && in2 && in3
        splittri1out!(csg, v1, v2, v3, triangles)
    elseif !in2 && in1 && in3
        splittri1out!(csg, v2, v3, v1, triangles)
    elseif !in3 && in1 && in2
        splittri1out!(csg, v3, v1, v2, triangles)
    end
end

function intersecttri!(csg::CSGTree{T}, tri::Triangle{T}, triangles::Vector{Triangle{T}}, thisforcoplanar = false) where {T<:Real}
    v1 = vertex(tri, 1)
    v2 = vertex(tri, 2)
    v3 = vertex(tri, 3)
    in1 = inside(csg, v1)
    in2 = inside(csg, v2)
    in3 = inside(csg, v3)
    if thisforcoplanar
        # if we are allowing dupes then we want to count the verts on the surface as being inside too
        # in case we have coplanar faces
        in1 |= onsurface(csg, v1)
        in2 |= onsurface(csg, v2)
        in3 |= onsurface(csg, v3)
    end
    if (in1 && in2 && in3)
        # whole triangle is valid
        push!(triangles, tri)
    elseif in1 && in2 && !in3
        splittri2out!(csg, v1, v2, v3, triangles)
    elseif in2 && in3 && !in1
        splittri2out!(csg, v2, v3, v1, triangles)
    elseif in1 && in3 && !in2
        splittri2out!(csg, v3, v1, v2, triangles)
    elseif in1 && !in2 && !in3
        splittri1out!(csg, v1, v2, v3, triangles)
    elseif in2 && !in1 && !in3
        splittri1out!(csg, v2, v3, v1, triangles)
    elseif in3 && !in1 && !in2
        splittri1out!(csg, v3, v1, v2, triangles)
    end
end

function splittri1out!(csg::CSGTree{T}, outv::SVector{3,T}, inv1::SVector{3,T}, inv2::SVector{3,T}, triangles::Vector{Triangle{T}}, recurse::Int = 0) where {T<:Real}
    if !validtri(outv, inv1, inv2) || recurse > VIS_RECURSION_LIMIT
        # if this triangle is invalid then just stop
        return
    end
    n1 = norm(outv - inv1)
    n2 = norm(outv - inv2)
    n3 = norm(inv1 - inv2)
    if min(n1, n2, n3) < MIN_VIS_TRI_SIZE
        # triangle is too small to bother splitting so just add it as is
        push!(triangles, Triangle(outv, inv1, inv2))
        return
    end
    n1 *= 1.0 + MESH_PRECISION
    n2 *= 1.0 + MESH_PRECISION
    int1 = closestintersection(evalcsg(csg, Ray(outv, inv1 - outv)), false)
    int2 = closestintersection(evalcsg(csg, Ray(outv, inv2 - outv)), false)
    if !(int1 === nothing || int2 === nothing) && validtri(outv, point(int1), point(int2)) && α(int1) <= n1 && α(int2) <= n2
        mid = (inv1 + inv2) / 2
        test = closestintersection(evalcsg(csg, Ray(outv, mid - outv)), false)
        if test !== nothing
            if !isapprox(point(test), mid, atol = MESH_PRECISION)
                splittri1out!(csg, outv, point(int1), point(test), triangles, recurse + 1)
                splittri1out!(csg, outv, point(test), point(int2), triangles, recurse + 1)
            else
                push!(triangles, Triangle(outv, point(int1), point(int2)))
            end
        else
            push!(triangles, Triangle(outv, point(int1), point(int2)))
        end
    end
end

function splittri2out!(csg::CSGTree{T}, outv1::SVector{3,T}, outv2::SVector{3,T}, inv::SVector{3,T}, triangles::Vector{Triangle{T}}, flip::Bool = false, recurse::Int = 0) where {T<:Real}
    if !validtri(outv1, outv2, inv) || recurse > VIS_RECURSION_LIMIT
        # if this triangle is invalid then just stop
        return
    end
    n1 = norm(outv1 - inv)
    n2 = norm(outv2 - inv)
    n3 = norm(outv1 - outv2)
    if min(n1, n2, n3) < MIN_VIS_TRI_SIZE
        # triangle is too small to bother splitting so just add it as is
        push!(triangles, Triangle(outv1, outv2, inv))
        return
    end
    n1 *= 1.0 + MESH_PRECISION
    n2 *= 1.0 + MESH_PRECISION
    int1 = closestintersection(evalcsg(csg, Ray(outv1, inv - outv1)), false)
    int2 = closestintersection(evalcsg(csg, Ray(outv2, inv - outv2)), false)
    if int2 !== nothing && α(int2) <= n2
        # int2 is valid
        if int1 !== nothing && α(int1) <= n1
            if isapprox(point(int1), point(int2), atol = MESH_PRECISION)
                # we just need one triangle in this case as inv lies on the boundary of the shape
                if validtri(outv1, outv2, point(int2))
                    if flip
                        push!(triangles, Triangle(outv1, point(int2), outv2))
                    else
                        push!(triangles, Triangle(outv1, outv2, point(int2)))
                    end
                end
            else
                if flip
                    splittri1out!(csg, outv1, point(int1), point(int2), triangles)
                else
                    splittri1out!(csg, outv1, point(int2), point(int1), triangles)
                end
                splittri2out!(csg, outv2, outv1, point(int2), triangles, !flip, recurse + 1)
            end
        else
            if validtri(outv1, outv2, point(int2))
                if flip
                    push!(triangles, Triangle(outv1, point(int2), outv2))
                else
                    push!(triangles, Triangle(outv1, outv2, point(int2)))
                end
            end
        end
    elseif int1 !== nothing && α(int1) <= n1
        if validtri(outv1, outv2, point(int1))
            if flip
                push!(triangles, Triangle(outv1, point(int1), outv2))
            else
                push!(triangles, Triangle(outv1, outv2, point(int1)))
            end
        end
    end
end

# visualization functions from CSG objects, most actualy work happens in the above functions
function makemeshi(a::UnionNode{T}, subdivisions::Int = 30)::Vector{Triangle{T}} where {T<:Real}
    # do the two children in parallel
    th1 = Threads.@spawn makemesh(a.leftchild, subdivisions)
    th2 = Threads.@spawn makemesh(a.rightchild, subdivisions)
    ltris = (fetch(th1)::TriangleMesh{T}).triangles
    rtris = (fetch(th2)::TriangleMesh{T}).triangles
    triangles = newintrianglepool!(T)
    @inbounds for tri in ltris
        uniontri!(a.rightchild, tri, triangles, true)
    end
    @inbounds for tri in rtris
        uniontri!(a.leftchild, tri, triangles)
    end
    return triangles
end

function makemeshi(a::IntersectionNode{T}, subdivisions::Int = 30)::Vector{Triangle{T}} where {T<:Real}
    # do the two children in parallel
    th1 = Threads.@spawn makemesh(a.leftchild, subdivisions)
    th2 = Threads.@spawn makemesh(a.rightchild, subdivisions)
    ltris = (fetch(th1)::TriangleMesh{T}).triangles
    rtris = (fetch(th2)::TriangleMesh{T}).triangles
    triangles = newintrianglepool!(T)
    @inbounds for tri in ltris
        intersecttri!(a.rightchild, tri, triangles, true)
    end
    @inbounds for tri in rtris
        intersecttri!(a.leftchild, tri, triangles)
    end
    return triangles
end

function makemeshi(a::ComplementNode{T}, subdivisions::Int = 30)::Vector{Triangle{T}} where {T<:Real}
    child_triangles = makemeshi(a.child, subdivisions)
    # flip the normals
    triangles = newintrianglepool!(T)
    @inbounds for i in 1:length(child_triangles)
        tri = child_triangles[i]
        push!(triangles, Triangle(vertex(tri, 3), vertex(tri, 2), vertex(tri, 1)))
    end
    return triangles
end

function makemeshi(a::LeafNode{T}, subdivisions::Int = 30)::Vector{Triangle{T}} where {T<:Real}
    child_triangles = triangulate(a.geometry, subdivisions)
    triangles = newintrianglepool!(T)
    @inbounds for i in 1:length(child_triangles)
        push!(triangles, a.transform * child_triangles[i])
    end
    return triangles
end

function makemesh(c::CSGTree{T}, subdivisions::Int = 30)::TriangleMesh{T} where {T<:Real}
    m = TriangleMesh(copy(makemeshi(c, subdivisions))) # need to copy with multithreading
    emptytrianglepool!(T)
    emptyintervalpool!(T)
    return m
end


text:
# --------------------------------------------------------------------------
# ACE.jl and SHIPs.jl: Julia implementation of the Atomic Cluster Expansion
# Copyright (c) 2019 Christoph Ortner <christophortner0@gmail.com>
# All rights reserved.
# --------------------------------------------------------------------------

using Documenter, ACE

makedocs(sitename="ACE.jl Documentation",
         pages = [
        "Home" => "index.md",
        "Introduction" => "intro.md",
        "Getting Started" => "gettingstarted.md",
        "Developer Docs" => "devel.md",
        "ED-Bonds" => "envpairbasis.md"
        # "Subsection" => [
        #     ...
        # ]
        ])

# deploydocs(
#     repo = "github.com/JuliaMolSim/ACE.jl.git",
# )


text:
# #  Copyright 2016 Eric S. Tellez <eric.tellez@infotec.mx>
# #
# #    Licensed under the Apache License, Version 2.0 (the "License");
# #    you may not use this file except in compliance with the License.
# #    You may obtain a copy of the License at
# #
# #        http://www.apache.org/licenses/LICENSE-2.0
# #
# #    Unless required by applicable law or agreed to in writing, software
# #    distributed under the License is distributed on an "AS IS" BASIS,
# #    WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# #    See the License for the specific language governing permissions and
# #    limitations under the License.

# export Is2014Search

# mutable struct Is2014Search <: LocalSearchAlgorithm
# end

# # From Malkov et al. 2014 Information System.
# # def KNNSearch(q, m, k) {
# #     TreeSet [object] tempRes, candidates, visitedSet, result
# #
# #     for (i = 0; i < m; i++) {
# #        put random entry point in candidates
# #        tempRes <-- null;
# #        repeat {
# #        	get element c closest from candidates to q;
# # 	        remove c from candidates;
# # 	        // check stop condition:
# # 	        if (c is further than k-th element from result) {
# # 	             break repeat;
# #     	    }
# # 	        // update list of candidates:
# # 	        for (every element e from friends of c) {
# # 	             if (e is not in visitedSet) {
# # 	                  add e to visitedSet, candidates, tempRes;
# # 		         }
# #           }
# # 	     }
# #    	 // aggregate the results:
# # 	     add objects from tempRes to result;
# #    }
# #    return best k elements from result;
# # }

# function greedy_search_is2014(index::SearchGraph{T}, q::T, res::R, tabu::BitVector, candidates::KnnResult) where {T,R <: Result}
#     nodeID::Int32 = rand(1:length(index.db))

#     # if !in(nodeID, tabu)
#     if !tabu[nodeID]
#         d::Float32 = convert(Float32, index.dist(index.db[nodeID], q))
#         # push!(tabu, nodeID)
#         tabu[nodeID] = true
#         push!(res, nodeID, d)
#         push!(candidates, nodeID, d)
#     end

#     while length(candidates) > 0
#         best = popfirst!(candidates)
#         # cov = last(res).dist
#         if best.dist > covrad(res)
#             break
#         end

#         @inbounds for childID in index.links[best.objID]
#             # if !in(childID, tabu)
#             if !tabu[childID]
#                 d = convert(Float32, index.dist(index.db[childID], q))
#                 # push!(tabu, childID)
#                 tabu[childID] = true
#                 push!(res, childID, d)
#                 push!(candidates, childID, d)
#             end
#         end
#     end
# end

# function search(algorithm::Is2014Search, index::SearchGraph{T}, q::T, res::R) where {T,R <: Result}
#     n = length(index.db)
#     tabu = falses(n)
#     restarts = min(abs(index.restarts), n)
#     candidates = KnnResult(128)

#     for i=1:restarts
#         greedy_search_is2014(index, q, res, tabu, candidates)
#     end

#     return res
# end


text:
# These functions are still WIP

function noise_model(freq_hz, fknee_hz, rms, alpha)
    rms * (1 + fknee_hz / freq_hz)^alpha
end

function baseline_spectrum(freq_hz, baseline_period_s)
    while false

    end
end


text:
# Raytracer.jl
# Raytracing for the generation of photorealistic images in Julia
# Copyright (c) 2021 Samuele Colombo, Paolo Galli

# Parser of SceneLang

"""
    ValueLoc

Stores a Julia value and the source location of its declaration.
"""
struct ValueLoc
    value::Any
    loc::SourceLocation
end

function Base.show(io::IO, valueloc::ValueLoc)
    print(io, valueloc.value)
    printstyled(io, " @ ", valueloc.loc, color = :light_black)
end

function Base.show(io::IO, ::MIME"text/plain", valueloc::ValueLoc)
    print(io, valueloc.value)
    printstyled(io, " @ ", valueloc.loc, color = :light_black)
end

"""
    IdTableKey

Alias for `Union{Type{<:TokenValue}, LiteralType}`. Used as the key type for [`IdTable`](@ref).
"""
const IdTableKey = Union{Type{<:TokenValue}, LiteralType}

"""
    IdTable

Alias to `Dict{IdTableKey, Dict{Symbol, ValueLoc}}`.

Dictionary with all the variables set in a SceneLang script.
Keys represent the equivalent [`TokenValue`](@ref) or [`LiteralType`](@ref) type to the Julia type stored in [`ValueLoc`](@ref).
"""
const IdTable = Dict{IdTableKey, Dict{Symbol, ValueLoc}}

function Base.show(io::IO, table::IdTable)
    show(io, MIME("text/plain"), table)
end

function Base.show(io::IO, ::MIME"text/plain", table::IdTable)
    for (key, subdict) ∈ pairs(table)
        printstyled(io, "#### ", key, " ####", color = :blue, bold = true)
        println()
        for (id, value) ∈ pairs(subdict)
            printstyled(io, id, color = :green)
            println(io, " = ", value)
        end
        println(io)
    end
end

"""
    RendererSettings

Struct containing a renderer type and a `NamedTuple` of the named arguments needed for its construction.

Since a [`Renderer`](@ref) type cannot be directly stored into a [`Scene`](@ref) due to it needing at least the [`World`](@ref) argument,
we can use this struct to store everything else, ready to be constructed.
"""
struct RendererSettings
    type::Type{<:Renderer}
    kwargs::NamedTuple
end

"""
    TracerSettings

Struct containing a `NamedTuple` of the named arguments needed to construct an [`ImageTracer`](@ref).

Since a [`ImageTracer`](@ref) type cannot be directly stored into a [`Scene`](@ref) due to it needing at least the [`Camera`](@ref)
and [`HdrImage`](@ref) arguments, we can use this struct to store everything else, ready to be constructed.
"""
struct TracerSettings
    kwargs::NamedTuple
end

"""
    CameraOrNot

Alias for `Union{Camera, Nothing}`.

See also: [`Scene`](@ref)
"""
const CameraOrNot = Union{Camera, Nothing}
"""
    RendererOrNot

Alias for `Union{RendererSettings, Nothing}`.

See also: [`Scene`](@ref)
"""
const RendererOrNot = Union{RendererSettings, Nothing}
"""
    ImageOrNot

Alias for `Union{HdrImage, Nothing}`.

See also: [`Scene`](@ref)
"""
const ImageOrNot = Union{HdrImage, Nothing}
"""
    TracerOrNot

Alias for `Union{TracerSettings, Nothing}`.

See also: [`Scene`](@ref)
"""
const TracerOrNot = Union{TracerSettings, Nothing}

"""
    Scene

A mutable struct containing all the significant elements of a renderable scene and all the declared variables of the SceneLang script.

# Fields

- `variables::`[`IdTable`](@ref): stores all the variables declared in the script
- `world::`[`World`](@ref): stores all the spawned [`Shape`](@ref)s
- `lights::`[`Lights`](@ref): stores all the spawned [`PointLight`](@ref)s
- `image::`[`ImageOrNot`](@ref): stores either the [`HdrImage`](@ref) to impress or `nothing`
- `camera::`[`CameraOrNot`](@ref): stores either the [`Camera`](@ref) or `nothing`
- `renderer::`[`RendererOrNot`](@ref): stores either the [`RendererSettings`](@ref) for the [`Renderer`](@ref) or `nothing`
- `tracer::`[`TracerOrNot`](@ref): strores either the [`TracerSettings`](@ref) for the [`ImageTracer`](@ref) or `nothing`
- `time::Float32`: stores the animation time
"""
mutable struct Scene
    variables::IdTable
    world::World
    lights::Lights
    image::ImageOrNot
	camera::CameraOrNot
	renderer::RendererOrNot
    tracer::TracerOrNot
    time::Float32
end

"""
    Scene(; variables::IdTable = IdTable(),
            world::World = World(),
            lights::Lights = Lights(),
            image::ImageOrNot = nothing,
            camera::CameraOrNot = nothing,
            renderer::RendererOrNot = nothing,
            tracer::TracerOrNot = nothing,
            time::Float32 = 0f0)

Constructor for a [`Scene`](@ref) instance.
"""
function Scene(; variables::IdTable = IdTable(),
                 world::World = World(),
                 lights::Lights = Lights(),
                 image::ImageOrNot = nothing,
                 camera::CameraOrNot = nothing,
                 renderer::RendererOrNot = nothing,
                 tracer::TracerOrNot = nothing,
                 time = 0f0)
	Scene(variables, world, lights, image, camera, renderer, tracer, time)
end

"""
    Scene(variables::Vector{Pair{Type{<:TokenValue}, Vector{Pair{Symbol, Token}}}};
          world::World = World(),
          lights::Lights = Lights(),
          image::ImageOrNot = nothing,
          camera::CameraOrNot = nothing,
          renderer::RendererOrNot = nothing,
          tracer::TracerOrNot = nothing,
          time = 0f0)

Constructor for a [`Scene`](@ref) instance.

Slightly more convenient than the default constructor to manually construct in a Julia code.
"""
function Scene(variables::Vector{Pair{Type{<:TokenValue}, Vector{Pair{Symbol, Token}}}};
               world::World = World(), lights::Lights = Lights(),
               image::ImageOrNot = nothing,
               camera::CameraOrNot = nothing,
               renderer::RendererOrNot = nothing,
               tracer::TracerOrNot = nothing,
               time = 0f0)

	variables =  Dict(zip(first.(variables), (Dict(last(pair)) for pair ∈ variables)))
	Scene(variables, world, lights, image, camera, renderer, tracer, time)
end

function Base.show(io::IO, ::MIME"text/plain", scene::Scene)
    printstyled(io, ".VARIABLES\n\n", color = :magenta, bold= true)
    show(io, scene.variables)
    println(io)

    printstyled(io, ".WORLD\n\n", color = :magenta, bold= true)
    if !isempty(scene.world)
        for (i, shape) ∈ enumerate(scene.world)
            printstyled(io, i, color = :green)
            println(io, ": ", shape)
        end
        println(io)
    end
    println(io)

    printstyled(io, ".LIGHTS\n\n", color = :magenta, bold= true)
    if !isempty(scene.lights)
        for (i, light) ∈ enumerate(scene.lights)
            printstyled(io, i, color = :green)
            println(io, ": ", light)
        end
        println(io)
    end
    println(io)

    printstyled(io, ".IMAGE\n\n", color = :magenta, bold= true)
    if !isnothing(scene.image)
        printstyled(io, "image", color = :green)
        println(io, " = ", scene.image)
        println(io)
    end
    println(io)

    printstyled(io, ".CAMERA\n\n", color = :magenta, bold= true)
    if !isnothing(scene.camera)
        printstyled(io, "camera", color = :green)
        println(io, " = ", scene.camera)
        println(io)
    end
    println(io)

    printstyled(io, ".RENDERER\n\n", color = :magenta, bold= true)
    if !isnothing(scene.renderer)
        printstyled(io, "type", color = :green)
        println(io, " = ", scene.renderer.type)
        if isempty(scene.renderer.kwargs)
            printstyled(io, "default kwargs\n\n", color=:light_black)
        else
            for (kw, value) ∈ pairs(scene.renderer.kwargs)
                printstyled(io, kw, color = :green)
                println(io, " = ", value)
            end
        end
    end

    printstyled(io, ".TRACER\n\n", color = :magenta, bold= true)
    if !isnothing(scene.tracer)
        if isempty(scene.tracer.kwargs)
            printstyled(io, "default kwargs\n\n", color=:light_black)
        else
            for (kw, value) ∈ pairs(scene.tracer.kwargs)
                printstyled(io, kw, color = :green)
                println(io, " = ", value)
            end
        end
    end

    printstyled(io, ".TIME\n\n", color = :magenta, bold= true)
    printstyled(io, "time", color = :green)
    println(io, " = ", scene.time)
    println(io)
end


text:
@testset "fft-tests" begin
    x = reshape(1:6, 2, 3)
    ax = AxisArray(x, 2:3, 3.0:5.0)

    @testset "fft" begin
        x_fft = fft(x, 2)
        ax_fft = fft(ax, 2)
        @test x_fft == ax_fft
        @test keys.(axes(ax_fft)) == keys.(axes(ax))
    end

    @testset "ifft" begin
        x_ifft = ifft(x, 2)
        ax_ifft = ifft(ax, 2)
        @test ax_ifft == x_ifft
        @test keys.(axes(ax_ifft)) == keys.(axes(ax))
    end

    @testset "bfft" begin
        x_bfft = bfft(x, 2)
        ax_bfft = bfft(ax, 2)
        @test ax_bfft == x_bfft
        @test keys.(axes(ax_bfft)) == keys.(axes(ax))
    end
end


text:
function generate_cartesian_orbit(obj::OMMFile)
    return jcall(obj, "generateCartesianOrbit", CartesianOrbit, ())
end

function generate_keplerian_orbit(obj::OMMFile)
    return jcall(obj, "generateKeplerianOrbit", KeplerianOrbit, ())
end

function generate_spacecraft_state(obj::OMMFile)
    return jcall(obj, "generateSpacecraftState", SpacecraftState, ())
end

function generate_tle(obj::OMMFile)
    return jcall(obj, "generateTLE", TLE, ())
end

function get_b_star(obj::OMMFile)
    return jcall(obj, "getBStar", jdouble, ())
end

function get_classification_type(obj::OMMFile)
    return jcall(obj, "getClassificationType", jchar, ())
end

function get_element_set_number(obj::OMMFile)
    return jcall(obj, "getElementSetNumber", JString, ())
end

function get_ephemeris_type(obj::OMMFile)
    return jcall(obj, "getEphemerisType", jint, ())
end

function get_mean_motion(obj::OMMFile)
    return jcall(obj, "getMeanMotion", jdouble, ())
end

function get_mean_motion_dot(obj::OMMFile)
    return jcall(obj, "getMeanMotionDot", jdouble, ())
end

function get_mean_motion_dot_dot(obj::OMMFile)
    return jcall(obj, "getMeanMotionDotDot", jdouble, ())
end

function get_meta_data(obj::OMMFile)
    return jcall(obj, "getMetaData", OMMFile_OMMMetaData, ())
end

function get_norad_id(obj::OMMFile)
    return jcall(obj, "getNoradID", Integer, ())
end

function get_rev_at_epoch(obj::OMMFile)
    return jcall(obj, "getRevAtEpoch", jint, ())
end

function get_tle_related_parameters_comment(obj::OMMFile)
    return jcall(obj, "getTLERelatedParametersComment", List, ())
end