text:
export
    AtiyahBottFormula

const  number_trees = [1, 1, 2, 3, 6, 11, 23, 47, 106, 235, 551, 1301, 3159]
#the number of non-isomorphic graphs with given number of vertices (starting from 2)

# the following function performs the computation of the Atiyah-Bott formula
# for a particular graph.

"""
    AtiyahBottFormula(n, d, m, P; do_check, show_bar, down_col)

Apply the Atiyah-Bott residue formula to the class `P`, in the moduli space of rational marked stable maps to the projective space of dimension `n` of degree `d` with `m` marks.
# Arguments
- `n::Int64`: the dimension of the projective space, it must be between ``1`` and ``254``.
- `d::Int64`: the degree of the stable maps, it must be between ``1`` and ``13``.
- `m::Int64`: the number of marks.
- `P`: the equivariant class.
- `do_check::Bool`: if `true`, checks if `P` is a well defined zero cycle, and stops the computation if this is not true. If `false`, the computation may have an unexpected behaviour. By default is `true`.
- `show_bar::Bool`: hide the progress bar if and only if this condition is `false`. By default is `true`.
- `down_col::Bool`: check if all colorations needed in the computation are in the folder Data, and download them otherwise. Once downloaded, those files can be used for future computations. By default is `true`.

The general construction of `P` is the following:
```julia-repl
julia> P =
```
After `=`, one has to write an expression in the equivariant classes. The expression is a combination of the equivariant classes. We compute the degree of `P` by

```julia-repl
julia> AtiyahBottFormula(n,d,m,P);
```

# Example
```julia-repl
julia> P = Hypersurface(5);
julia> AtiyahBottFormula(3,1,0,P);
Warning: the class is not a 0-cycle.
julia> AtiyahBottFormula(4,1,0,P);
Result: 2875
julia> AtiyahBottFormula(4,1,0,P,do_check = false); #skip the preliminary check on `P`
julia> AtiyahBottFormula(4,1,0,P,show_bar = false); #it does not show the progress bar
julia> AtiyahBottFormula(4,1,0,P,down_col = false); #it does not download the colorations, if missing
```

The function returns an array of the same dimension of `P` (non-vectorized classes are assumed as 1-dimensional arrays). The Julia notation for accessing to array is `name_of_array[i]` where `i` is an index starting from 1.

# Example
```julia-repl
julia> P = Incidency(2)*Hypersurface(3);
julia> x = AtiyahBottFormula(3,2,0,P)[1];
Result: 81
julia> x
81
```

The class `P` supports parameters.
```julia-repl
julia> P = Hypersurface(3)*(Incidency(2)//3)^(d-1);
julia> d = 2;
julia> AtiyahBottFormula(3,d,0,P);
Result: 27
julia> d = 3;
julia> AtiyahBottFormula(3,d,0,P);
Result: 84
```

More examples are available in the support of the equivariant classes. It is enough to type `?` and then the name of the class. Currently, the supported classes are:

* `O1_i`         (Euler class of ``\\mathrm{ev}_i^*\\mathcal{O}_{\\mathbb{P}^n}(1)``)
* `O1`           (product of all `O1_i`)
* `Psi`          (cycle of ``\\psi``-classes)
* `Jet`          (Euler class of the jet bundle ``J^p`` of ``\\mathrm{ev}^*\\mathcal{O}_{\\mathbb{P}^n}(q)``)
* `Hypersurface` (Euler class of the direct image of ``\\mathrm{ev}^*\\mathcal{O}_{\\mathbb{P}^n}(b)``)
* `Incidency`    (cycle parameterizing curves meeting a linear subspace)
* `Contact`      (cycle parameterizing contact curves)
* `R1`           (first derived functor of direct image of ``\\mathrm{ev}^*\\mathcal{O}_{\\mathbb{P}^n}(-k)``)

To add more classes, please contact the authors.
"""
function AtiyahBottFormula(n::Int64, deg::Int64, n_marks::Int64, P_input; do_check::Bool = true, show_bar::Bool = true, down_col::Bool = true)::Vector{fmpq}

    if n < 1
        printstyled("ERROR: ", bold=true, color=:red)
        println("n must be positive, correct ", n)
        return [fmpq(0)]
    end
    if deg > 13 || deg < 1
        printstyled("ERROR: ", bold=true, color=:red)
        println("d must be between 1 and 13, correct ", deg)
        return [fmpq(0)]
    end
    if n_marks < 0
        printstyled("ERROR: ", bold=true, color=:red)
        println("m must be non negative, correct ", n_marks)
        return [fmpq(0)]
    end

    local n_results::Int64 = 0

    if isa(P_input, Array)
        n_results = length(P_input)
    else
        n_results = 1
    end
    local P::Vector{Function} = Vector(undef, n_results)

    if isa(P_input, Array)
        for i in eachindex(P)
            P[i] = P_input[i].func
        end
    else
        P[1] = P_input.func
    end

    if do_check && !is_zero_cycle(n, deg, n_marks, P)
        return [fmpq(0)]
    end

    if down_col && !fill_Data(n, deg)
        printstyled("ERROR: ", bold=true, color=:red)
        println("I was unable to download the colorations from the repository. Check your internet connection or execute:")
        printstyled("julia> ", bold=true, color=:light_green)
        println("AtiyahBottFormula($n,$deg,$n_marks,P,down_col = false);")
        return [fmpq(0)]
    end

    # local n_results::Int64 = length(P)   #this store the number of final results of our computation
    local result::Vector{Vector{fmpq}} = [[fmpq() for _ in 1:n_results] for _ in 1:Threads.nthreads()]
    local s::Vector{fmpq} = fmpq.(rand(UInt32, n+1))

    if show_bar #set up progress data
        threshold::Int64 = sum(v -> number_trees[v-1]*(n+1)*(n^(v-1))*(v^n_marks), 2:deg+1)
        progress_bar::Progress = Progress(threshold, barglyphs=BarGlyphs("[=> ]"), color=:green)
        current_graph::Threads.Atomic{Int64} = Threads.Atomic{Int}(0)
    end

    #Create the string_array

    local all_weights::Vector{Vector{Vector{Int}}} = Vector{Vector{Vector{Int}}}(undef, deg)

    list_g::IOStream = open(current_dir*"/list_trees.txt", "r")
    #open the file containing the list of Prufer sequences of graphs

    for ne in 1:deg #run the computation among all graphs with fixed number of vertices

        all_weights[ne] = get_weights(ne, deg)

    end

    dir::String = dirname(current_dir)*"/Data/Dimension$n/"
    read_dir::Union{Vector{String},Nothing} = nothing

    if !down_col
        if isdir(dir)
            read_dir = readdir(dirname(current_dir)*"/Data/Dimension$n/")
        else
            read_dir = String[]
        end
    end

    graph_generator = graph_gen(list_g, sum(number_trees[1:deg]), read_dir)

    l = Threads.SpinLock()

    Threads.@threads for str in collect(graph_generator)

        pruf_str = str[1]
        top_aut = str[2]
        local g::SimpleGraph = PruferToGraph(UInt8[Read[string(pruf_str[i])] for i in 1:length(pruf_str)])

        local c::Vector{UInt8} = Vector{UInt8}(undef, nv(g))
        local aut::Int64 = top_aut

        if str[3]
            file_name = dir*pruf_str*"0.gz"
            local cols = graph_coloring_from_file(GzipDecompressorStream(open(file_name, "r")))
        else
            local cols = graph_coloring(g, UInt8(n+1))
        end

        for col in cols   #we run among all colorations of g

            if str[3]
                c, aut = col
            else
                c = col
            end

            for m in Marks(nv(g), n_marks)    #we run among all marks of g, if n_marks==0 we have only the empty mark

                for w in all_weights[ne(g)]          #we run among all weights of g
                    try
                        local Euler::fmpq = fmpq(0)

                        eq!(Euler, Euler_inv(g,c,w,s,m))
                        div_eq!(Euler, aut*prod(w))

                        for res in 1:n_results      #compute each term of the array P
                            local temp::fmpq = fmpq(0)
                            eq!(temp, P[res](g,c,w,s,m))
                            mul_eq!(temp, Euler)
                            add_eq!(result[Threads.threadid()][res], temp)
                            # result[res] += P[res](g,c,w,s,m)*Euler    #apply Atiyah-Bott
                        end

                    catch err
                        if isa(err, DivideError)
                            error("Some division by zero occurred. Try again")
                        end
                        println(err)
                        error("Some error occurred")
                        return zeros(fmpq, n_results)
                    end
                end

                if show_bar
                    Threads.atomic_add!(current_graph, top_aut÷aut)
                    #progress_data.current_graph += progress_data.top_aut÷aut
                    #update the progress bar
                    Threads.lock(l)
                    update!(progress_bar, current_graph[],
                            showvalues = [(:"Total number of graphs",threshold),(:"Current graph",current_graph[])])
                    Threads.unlock(l)
                end
            end
        end
    end

    for nt in 2:Threads.nthreads()
        add_eq!(result[1], result[nt])
    end

    if n_results == 1
        println("Result: ", result[1][1])
    else
        for res in 1:n_results
            println("Result number ", res, ": ", result[1][res])
        end
    end

    return result[1]
end


text:
module VarianceComponentTypeTest

using VarianceComponentModels, MathProgBase, Ipopt, LinearAlgebra, Random, Test

Random.seed!(123)

# generate data from a d-variate response variane component model
n = 100   # no. observations
d = 2     # no. categories
m = 2     # no. variance components
p = 2     # no. covariates
X = randn(n, p)
B = ones(p, d)
Σ = ntuple(x -> zeros(d, d), m)
for i in 1:m
  Σi = randn(d, d)
  copyto!(Σ[i], Σi' * Σi)
end
## make the first variance component 0 matrix
#fill!(Σ[1], 0.0)
V = ntuple(x -> zeros(n, n), m)
for i = 1:m-1
  Vi = randn(n, 50)
  copyto!(V[i], Vi * Vi')
end
copyto!(V[m], Matrix(I, n, n))
# form Ω
Ω = zeros(n*d, n*d)
for i = 1:m
  global Ω += kron(Σ[i], V[i])
end
Ωchol = cholesky(Ω)
Y = X * B + reshape(Ωchol.L * randn(n*d), n, d)

@info "VarianceComponentModel type"
#@code_warntype VarianceComponentModel(B, Σ)
@inferred VarianceComponentModel(B, Σ)
#@code_warntype VarianceComponentModel(Σ)
@inferred VarianceComponentModel(Σ)
vcmodel = VarianceComponentModel(Σ)
@test isempty(vcmodel.B)

@info "TwoVarCompModelRotate type"
F = eigen(Σ[1], Σ[2])
#@code_warntype TwoVarCompModelRotate(B * F.vectors, F.values, F.vectors, logdet(Σ[2]))
@inferred TwoVarCompModelRotate(B * F.vectors, F.values, F.vectors, logdet(Σ[2]))

@info "VarianceComponentVariate type"
#@code_warntype VarianceComponentVariate(Y, V)
@inferred VarianceComponentVariate(Y, V)
vcdata = VarianceComponentVariate(Y, V)
@test isempty(vcdata.X)

@info "Forming VarianceComponentModel from data"
vcdata = VarianceComponentVariate(Y, X, V)
#@code_warntype VarianceComponentVariate(Y, X, V)
@inferred VarianceComponentVariate(Y, X, V)
vcdata = VarianceComponentVariate(Y, X, V)
#@code_warntype VarianceComponentModel(vcdata)
@inferred VarianceComponentModel(vcdata)
vcmodel = VarianceComponentModel(vcdata)
@test norm(vcmodel.B) ≈ 0.0
@test norm(vcmodel.Σ[1] - Matrix(I, d, d)) ≈ 0.0
@test norm(vcmodel.Σ[2] - Matrix(I, d, d)) ≈ 0.0

@info "Pre-compute eigen-decomposition and rotate data"
#@code_warntype TwoVarCompVariateRotate(vcdata)
@inferred TwoVarCompVariateRotate(vcdata)
vcdatarot = TwoVarCompVariateRotate(vcdata)
#@code_warntype TwoVarCompModelRotate(vcmodel)
@inferred TwoVarCompModelRotate(vcmodel)
vcmodelrot = TwoVarCompModelRotate(vcmodel)

@info "Test generalized eigen-decomposition"
vcdatacopy = deepcopy(vcdata)
rmul!(vcdatacopy.V[2], 2.0) # V[2] = 2I
vcdatacopyrot = TwoVarCompVariateRotate(vcdatacopy)
@test vcdatacopyrot.logdetV2 ≈ n * log(2.0) # det(V[2]) = 2^n

@info "VarianceComponentModel from TwoVarCompVariateRotate"
@inferred VarianceComponentModel(vcdatarot)
vcmodelfromrot = VarianceComponentModel(vcdatarot)
@test size(vcmodelfromrot.B, 1) == p
@test size(vcmodelfromrot.B, 2) == d
@test length(vcmodelfromrot.Σ) == m

@info "VarianceComponentAuxData from VarianceComponentVariate"
vcobsaux1 = VarianceComponentAuxData(vcdata)
vcobsaux2 = VarianceComponentAuxData(vcdatarot)
@test size(vcobsaux1.res) == size(vcobsaux2.res)
@test size(vcobsaux1.Xwork) == size(vcobsaux2.Xwork)
@test size(vcobsaux1.ywork) == size(vcobsaux2.ywork)
@test size(vcobsaux1.obswt) == size(vcobsaux2.obswt)

@info "Query functions"
@test eltype(vcmodel) == eltype(B)
@test eltype(vcdata) == eltype(Y)
@test eltype(vcmodelrot) == eltype(vcmodel)
@test eltype(vcdatarot) == eltype(vcdata)
@test length(vcmodel) == d
@test length(vcdata) == d
@test length(vcmodelrot) == d
@test length(vcdatarot) == d
@test size(vcdata) == (n, d)
@test size(vcdatarot) == (n, d)
@test nvarcomps(vcmodel) == m
@test nvarcomps(vcdata) == m
@test nvarcomps(vcmodelrot) == m
@test nvarcomps(vcdatarot) == m
@test nmeanparams(vcmodel) == p * d
@test nmeanparams(vcmodelrot) == p * d
@test nmeanparams(vcdata) == p * d
@test nmeanparams(vcdatarot) == p * d
@test nvarparams(vcmodel) == m * binomial(d + 1, 2)
@test nvarparams(vcmodelrot) == m * binomial(d + 1, 2)
@test nparams(vcmodel) == p * d + m * binomial(d + 1, 2)
@test nparams(vcmodelrot) == p * d + m * binomial(d + 1, 2)

@info "Mean, covariance, and residual of model"
#@code_warntype mean(vcmodel, vcdata)
@inferred mean(vcmodel, vcdata)
μ = mean(vcmodel, vcdata)
@test size(μ, 1) == n
@test size(μ, 2) == d
#@code_warntype cov(vcmodel, vcdata)
@inferred cov(vcmodel, vcdata)
Ω = cov(vcmodel, vcdata)
@test size(Ω, 1) == n * d
@test size(Ω, 2) == n * d
#@code_warntype residual(vcmodel, vcdata)
@inferred residual(vcmodel, vcdata)
res = residual(vcmodel, vcdata)
@test size(res, 1) == n
@test size(res, 2) == d
@inferred residual(vcmodelrot, vcdatarot)
resrot = residual(vcmodelrot, vcdatarot)
@test size(resrot, 1) == n
@test size(resrot, 2) == d

end # module VarianceComponentTypeTest

ext: jl

max_stars_repo_name: Hua-Zhou/VarianceComponentModels

max_issues_repo_name: Hua-Zhou/VarianceComponentModels

max_forks_repo_name: Hua-Zhou/VarianceComponentModels


text:
#HW2 Problem 4
#Minhao Yan my497, Wentong Chen wc422

#1. Recast the problem as an unconstrained utility maximization problem
#Plug the budget constraint back to the utility function
#u(c1) = (κ*c1^((η-1)/η)+(1-κ)*((w-p1*c1)/p2)^((η-1)/η))^(η/(η-1))

#2. Solve this utility maximization problem
Pkg.add("Optim")
using Optim

function utility_maximizer(κ, η, p1, p2, w, initial_guess, solver)
    u(c1)=-((κ*c1^((η-1)/η)+(1-κ)*((w-p1*c1)/p2)^((η-1)/η)))^(η/(η-1))
    results = optimize(u,initial_guess)
    c1_min=Optim.minimizer(results)[1]
    c2_min= (w-p1*c1_min)/p2
    u_max=-Optim.minimum(results)
    println("c1= $c1_min")
    println("c2= $c2_min")
    println("u_max= $u_max")
    return [c1_min c2_min u_max]
    return results
end

#3.Solve analytically
κ=.5
η=3
w=50
p1=2
p2=3
initial_guess=0.0
Result_ana = utility_maximizer(κ, η, p1, p2, w, initial_guess, Newton())

#Solve numerically
u1(c1)=κ*c1^((η-1)/η)+(1-κ)*((w-p1*c1)/p2)^((η-1)/η)
function g!(storage, x)
    storage[1]=-(η/(η-1))*u1(x[1])^(1/(η-1))*(κ*((η-1)/η)*x[1]^(-1/η)+(1-κ)*((η-1)/η)*((w-p1*x[1])/p2)^(-1/η)*(-p1/p2))
    storage[2] = 0.0
end
function utility_maximizer_num(κ, η, p1, p2, w, initial_guess, solver)
    x=initial_guess
    u(c1)=0.0*x[2]-((κ*c1^((η-1)/η)+(1-κ)*((w-p1*c1)/p2)^((η-1)/η)))^(η/(η-1))
    results = optimize(u, g!, initial_guess)
    c1_min=Optim.minimizer(results)[1]
    c2_min= (w-p1*c1_min)/p2
    u_max=-Optim.minimum(results)
    println("c1= $c1_min")
    println("c2= $c2_min")
    println("u_max= $u_max")
    return [c1_min c2_min u_max]
    return results
end

κ=.5
η=3
w=50
p1=2
p2=3
initial_guess=0.0
Return_num = utility_maximizer_num(κ, η, p1, p2, w, initial_guess, Newton())


text:
##############################################################################
#
# Rayleigh distribution from Distributions Handbook
#
##############################################################################

immutable Rayleigh <: ContinuousUnivariateDistribution
    scale::Float64
    function Rayleigh(s::Real)
        if s > 0.0
            new(float64(s))
        else
            error("scale must be positive")
        end
    end
end

Rayleigh() = Rayleigh(1.0)

insupport(d::Rayleigh, x::Number) = isreal(x) && isfinite(x) && 0.0 < x

kurtosis(d::Rayleigh) = d.scale^4 * (8.0 - ((3.0 * pi^2) / 4.0))

mean(d::Rayleigh) = d.scale * sqrt(pi / 2.)

median(d::Rayleigh) = d.scale * sqrt(2. * log(2.))

function pdf(d::Rayleigh, x::Real)
    if insupport(d, x)
        return (x / (d.scale^2)) * exp(-(x^2)/(2.0 * (d.scale^2)))
    else
        return 0.0
    end
end

rand(d::Rayleigh) = d.scale * sqrt(-2.0 * log(rand()))

skewness(d::Rayleigh) = d.scale^3 * (pi - 3.0) * sqrt(pi / 2.0)

var(d::Rayleigh) = d.scale^2 * (2.0 - pi / 2.0)


text:
using BenchmarkTools
using NumCME

function explore_states(state_type::Type{ST}, index_type::Type{T}) where {ST<:Integer, T<: Integer}
    S = Matrix{state_type}([[1;0] [-1;0] [0;1] [0;-1]])
    statespace = StateSpaceSparse(S, Vector{state_type}([0,0]); index_type=index_type)
    expand!(statespace, 200)
    return statespace
end

@btime space16 = explore_states(Int32, UInt32);
@btime space64 = explore_states(Int64, UInt64);


text:
#=
Functionality regarding spatial coordinates of an array, i.e. Longitude and Lattitude
as well as equal area grid handling and anything else related to space.
=#
using Statistics, StatsBase, StaticArrays
export SVector # for equal area grid
#########################################################################
# Spatial indexing
#########################################################################
export spatialidxs, lonlatfirst

"""
    spatialidxs(A::ClimArray) → idxs
Return an iterable that can be used to access all spatial points of `A` with the syntax
```julia
idxs = spatialidxs(A)
for i in idxs
    slice_at_give_space_point = A[i...]
end
```
Works for all types of space (`...` is necessary because `i` is a `Tuple`).
"""
spatialidxs(A::AbDimArray) = spatialidxs(spacestructure(A), A)
function spatialidxs(::OrthogonalSpace, A)
    lons = (Lon(i) for i in 1:size(A, Lon))
    lats = (Lat(i) for i in 1:size(A, Lat))
    return Iterators.product(lons, lats)
end

"""
    lonlatfirst(A::ClimArray, args...) → B
Permute the dimensions of `A` to make a new array `B` that has first dimension longitude,
second dimension latitude, with the remaining dimensions of `A` following
(useful for most plotting functions). Optional extra dimensions
can be given as `args...`, specifying a specific order for the remaining dimensions.

Example:
```julia
B = lonlatfirst(A)
C = lonlatfirst(A, Time)
```
"""
function lonlatfirst(C, args...)
    permutedims(C, (Lon, Lat, args..., otherdims(C, (Lon, Lat, args...))...))
end

#########################################################################
# Periodicity of longitude
#########################################################################
export lon_distance, wrap_lon, longitude_circshift
"""
    wrap_lon(x)
Wrap given longitude to -180 to 180 degrees.
"""
wrap_lon(x) = @. -180 + (360 + ((x+180) % 360)) % 360

"""
    lon_distance(λ1, λ2, Δλ = 360) → δ
Calculate distance `δ` (also in degrees) between longitudes `λ1, λ2`, but taking into
account the periodic nature of longitude, which has period `Δλ = 360ᵒ`.
"""
function lon_distance(x, y, p = eltype(x)(360))
    moddis = mod(abs(x - y), p)
    min(moddis, p - moddis)
end

"""
    longitude_circshift(X::ClimArray [, l]; wrap = true) → Y::ClimArray
Perform the same action as `Base.circshift`, but only for the longitudinal dimension
of `X` with shift `l`. If `wrap = true` the longitudes are wrapped to (-180, 180) degrees
using the modulo operation.

If `l` is not given, it is as much as necessary so that all longitudes > 180 are
wrapped.
"""
function longitude_circshift(X::ClimArray, l = nothing; wrap = true)
    if isnothing(l); l = count(≥(180), dims(X, Lon).val); end
    isnothing(l) && return X
    shifts = map(d -> d isa Lon ? l : 0, dims(X))
    shifted_data = circshift(X.data, shifts)
    shifted_lon = circshift(dims(X, Lon).val, l)
    if wrap; shifted_lon = wrap_lon.(shifted_lon); end
    shifted_lon = vector2range(shifted_lon)
    shifted_dim = Lon(shifted_lon; metadata = DimensionalData.metadata(dims(X, Lon)))
    new_dims = map(d -> d isa Lon ? shifted_dim : d, dims(X))
    return ClimArray(shifted_data, new_dims; name = X.name, attrib = X.attrib)
end

#########################################################################
# averaging functions over space or time
#########################################################################
using StatsBase

export latmean, spacemean, zonalmean, spaceagg, uniquelats

"""
    zonalmean(A::ClimArray [, W])
Return the zonal mean of `A`. Works for both [`OrthogonalSpace`](@ref) as well as
[`CoordinateSpace`](@ref). Optionally provide statistical weights `W`.
These can be the same `size` as `A` or only having the same latitude structure as `A`.
"""
zonalmean(A::AbDimArray, W = nothing) = zonalmean(spacestructure(A), A, W)
zonalmean(::OrthogonalSpace, A::AbDimArray, W) = dropagg(mean, A, Lon, W)

"""
    latmean(A::ClimArray)
Return the latitude-mean `A` (mean across dimension `Lat`).
This function properly weights by the cosine of the latitude.
"""
function latmean(A::AbDimArray)
    lw = _latweights(A)
    if ndims(A) > 1
        return dropagg(sum, broadcast_dims(*, A, lw), Lat)
    else
        return sum(A .* lw)
    end
end

# Warning!!! `_latweights` divides by the weight sum, because it is intended to be
# used only with the `sum` function (for A)
_latweights(A::AbDimArray) = _latweights(dims(A, Lat))
function _latweights(A::Lat)
    we = cosd.(A.val)
    we ./= sum(we)
    return ClimArray(we, (A,))
end

using StatsBase

"""
    spacemean(A::ClimArray [, W]) = spaceagg(mean, A, W)
Average `A` over its spatial coordinates.
Optionally provide statistical weights in `W`.
"""
spacemean(A, exw=nothing) = spaceagg(mean, A, exw)

"""
    spaceagg(f, A::ClimArray, W = nothing)
Aggregate `A` using function `f` (e.g. `mean, std`) over all available space (i.e.
longitude and latitude) of `A`, weighting every part of `A` by its spatial area.

`W` can be extra weights, to weight each spatial point with. `W` can either be a
`ClimArray` with same spatial information as `A`, or having exactly same dimensions as `A`.
"""
spaceagg(f, A::AbDimArray, exw=nothing) = spaceagg(spacestructure(A), f, A, exw)
function spaceagg(::OrthogonalSpace, f, A::AbDimArray, w=nothing)
    wtype = spaceweightassert(A, w)
    cosweights = repeat(cosd.(dims(A, Lat).val)', size(A, Lon))
    if dimindex(A, Lon) > dimindex(A, Lat)
        error("At the moment this function assumes that Lon preceeds Lat, use `permutedims`.")
    end
    other = otherdims(A, (Lon, Lat))
    # pre-calculate weights if possible
    if wtype == :no
        W = weights(cosweights)
    elseif wtype == :d2
        W = weights(cosweights .* Array(w))
    end
    # if A is only 2-dimensional, operation is trivial and single number
    if ndims(A) == 2
        return f(A, W)
    end
    # do the weighted average
    oidxs = otheridxs(A, (Lon(), Lat()))
    if wtype != :dany
        r = map(i -> f(view(A, i), W), oidxs)
    else
        r = map(i -> f(view(A, i), weights(view(w, i) .* cosweights)), oidxs)
    end
    return ClimArray(r, other, A.name)
end

function spaceweightassert(A, w)
    if !isnothing(w)
        wdims = dims(w)
        # Check that longitude/latitude match even if the actual dimensions are not
        # identical (because they are e.g. not loaded from the same file)
        if length(wdims) == 2
            @assert Lat ∈ basetypeof.(wdims)
            @assert Lon ∈ basetypeof.(wdims)
            @assert val.(wdims) == val.(dims(A, (Lon, Lat)))
            wtype = :d2
        else
            @assert basetypeof.(wdims) == basetypeof.(dims(A))
            @assert size(w) == size(A)
            wtype = :dany
        end
    else
        wtype = :no
    end
    return wtype
end

#########################################################################
# Hemispheric
#########################################################################
export hemispheric_means, hemispheric_functions

"""
    hemispheric_functions(A::ClimArray) → north, south
Return two arrays `north, south`, by splitting `A` to its northern and southern hemispheres,
appropriately translating the latitudes of `south` so that both arrays have the same
latitudinal dimension (and thus can be compared and do opperations between them).
"""
hemispheric_functions(A) = hemispheric_functions(spacestructure(A), A)
function hemispheric_functions(::OrthogonalSpace, A)
    nh = A[Lat(0..90)]
    sh = A[Lat(-90..0)]
    # TODO: this can be a function "reverse dim"
    di = dimindex(sh, Lat)
    newdims = [dims(sh)...]
    newdims[di] = dims(nh, Lat)
    data = reverse(Array(sh); dims = di)
    sh = ClimArray(data, (newdims...,))
    return nh, sh
end

"""
    hemispheric_means(A [,W]) → nh, sh
Return the (proper) averages of `A` over the northern and southern hemispheres.
Notice that this function explicitly does both zonal as well as meridional averaging.
Use [`hemispheric_functions`](@ref) to just split `A` into two hemispheres.

Optionally provide weights `W` that need to have the same structure as [`spaceagg`](@ref).
"""
hemispheric_means(A, args...) = hemispheric_means(spacestructure(A), A, args...)
function hemispheric_means(::OrthogonalSpace, A::AbDimArray)
    @assert hasdim(A, Lat)
    if hasdim(A, Lon)
        B = zonalmean(A)
    else
        B = A
    end
    nh = latmean(B[Lat(0..90)])
    sh = latmean(B[Lat(-90..0)])
    return nh, sh
end

latitudes(A) = latitudes(spacestructure(A), A)
latitudes(::OrthogonalSpace, A) = dims(A, Lat).val
#########################################################################
# Tropics/extratropics
#########################################################################
export tropics_extratropics
"""
    tropics_extratropics(A::ClimArray; lower_lat=30, higher_lat=90) → tropics, extratropics
Separate the given array into two arrays: one having latitudes ℓ ∈ [-lower_lat, +lower_lat], and one
having [-higher_lat:-lower_lat, lower_lat:higher_lat].
"""
tropics_extratropics(A, args...; kwargs...) =
tropics_extratropics(spacestructure(A), A, args...; kwargs...)

function tropics_extratropics(::OrthogonalSpace, A; lower_lat=30, higher_lat=90)
    tropics = A[Lat(-lower_lat..lower_lat)]
    latdim = dims(A, Lat)
    extra_idxs_sh = DimensionalData.selectindices(latdim, -higher_lat..(-lower_lat))
    extra_idxs_nh = DimensionalData.selectindices(latdim, lower_lat..higher_lat)
    extra_idxs = vcat(extra_idxs_sh, extra_idxs_nh)
    extratropics = A[Lat(extra_idxs)]
    return tropics, extratropics
end


text:
# Section 4.4 of
# A. A. Ahmadi, and A. Majumdar
# DSOS and SDSOS Optimization: More Tractable Alternatives to Sum of Squares and Semidefinite Optimization
# 2017

using MultivariateMoments

@testset "Options Pricing with $(factory.constructor)" for factory in sdp_factories
    isscs(factory) && continue
    continue # see https://github.com/JuliaOpt/SumOfSquares.jl/issues/48
    @polyvar x y z
    σ = [184.04, 164.88, 164.88, 184.04, 164.88, 184.04]
    X = [x^2, x*y, x*z, y^2, y*z, z^2, x, y, z, 1]
    μ = measure([σ .+ 44.21^2; 44.21 * ones(3); 1],
                X)
    function optionspricing(K, cone)
        m = SOSModel(factory)
        @variable m p Poly(X)
        @constraint m p in cone
        @constraint m p - (x - K) in cone
        @constraint m p - (y - K) in cone
        @constraint m p - (z - K) in cone
        @objective m Min dot(μ, p)
        JuMP.optimize!(m)
        @test JuMP.primal_status(m) == MOI.FeasiblePoint
        JuMP.objective_value(m)
    end

    K = [30, 35, 40, 45, 50]
    dsos_codsos_exp   = [132.63, 132.63, 132.63, 132.63, 132.63]
    sdsos_cosdsos_exp = [ 21.51,  17.17,  13.20,   9.85,   7.30]
    sos_cosos_exp = sdsos_cosdsos_exp
    for i in eachindex(K)
        @test optionspricing(K[i], CoDSOSCone())  ≈ dsos_codsos_exp[i]   atol=1e-2
        @test optionspricing(K[i], CoSDSOSCone()) ≈ sdsos_cosdsos_exp[i] atol=1e-2
        @test optionspricing(K[i], CoSOSCone())   ≈ sos_cosos_exp[i]     atol=1e-2
    end
end


text:
# Activate environment

using ClusterManagers,Distributed

# Add in the cores allocated by the scheduler as workers
addprocs(SlurmManager(parse(Int,ENV["SLURM_NTASKS"])-1))
print("Added workers: ")
println(nworkers())

# Load required packages and helpers on all processes
@everywhere begin
    using DistributedArrays,Statistics
    include("../parallelhelpers.jl")
end

# Load the file names
dataLoc = "../../../../data/word_count/";
fnames = dataLoc.*readdir(dataLoc)

# Create a distributed array with the file names
dfnames = distribute(fnames)

# Calculate the word counts
# Broadcast over a distributed array automatically does the operation on the correct worker
allcounts =  getcounts.(dfnames)

# Collect the distributed array onto the leader process
allcountsL = convert(Array,allcounts)

# Calculate the overall counts on the leader
overallcounts = getoverallcount(allcountsL)

# Calculate the normalized counts on the leader and print results
normcount = getnormcount.(allcountsL)

max_stars_repo_name: llsc-supercloud/teaching-examples

max_issues_repo_name: llsc-supercloud/teaching-examples

max_forks_repo_name: llsc-supercloud/teaching-examples


text:
using Test
using NCDatasets
using CalibrateEDMF.KalmanProcessUtils
import CalibrateEDMF.KalmanProcessUtils: LearningRateScheduler

@testset "KalmanProcessUtils" begin

    @testset "LearningRateScheduler" begin
        Δt_init = 1.0
        τ = 5

        Δt_decay = PiecewiseConstantDecay(Δt_init, τ)
        Δt_growth = PiecewiseConstantGrowth(Δt_init, τ)
        @test isa(Δt_decay, LearningRateScheduler)
        @test isa(Δt_growth, LearningRateScheduler)

        t_future = 10
        @test get_Δt(Δt_decay, t_future) < Δt_init
        @test get_Δt(Δt_growth, t_future) > Δt_init
    end

end


text:
using WindMixing: plot_LES_3D
using Images, FileIO
using ImageTransformations
using CairoMakie

u_img = rotr90(load(assetpath(joinpath(pwd(), "images_plot\\u.png"))))
v_img = rotr90(load(assetpath(joinpath(pwd(), "images_plot\\v.png"))))
T_img = rotr90(load(assetpath(joinpath(pwd(), "images_plot\\T.png"))))
uw_img = rotr90(load(assetpath(joinpath(pwd(), "images_plot\\uw.png"))))
vw_img = rotr90(load(assetpath(joinpath(pwd(), "images_plot\\vw.png"))))
wT_img = rotr90(load(assetpath(joinpath(pwd(), "images_plot\\wT.png"))))
Ri_img = rotr90(load(assetpath(joinpath(pwd(), "images_plot\\Ri.png"))))
z_img = load(assetpath(joinpath(pwd(), "images_plot\\z.png")))
T_3D_img = rotr90(load(assetpath(joinpath(pwd(), "images_plot\\T_3D.png"))))

axis_images = (
    u = u_img,
    v = v_img,
    T = T_img,
    z = z_img,
    uw = uw_img,
    vw = vw_img,
    wT = wT_img,
    Ri = Ri_img,
    T_3D = T_3D_img,
)

# FILE_DIR = "C:\\Users\\xinle\\Downloads\\three_layer_constant_fluxes_linear_hr144_Qu5.5e-04_Qb3.0e-08_f1.0e-04_Nh256_Nz128_strong_wind_weak_cooling"
FILE_DIR = "C:\\Users\\xinle\\Downloads\\three_layer_constant_fluxes_linear_hr144_Qu0.0e+00_Qb5.0e-08_f1.0e-04_Nh256_Nz128_free_convection"
# FILE_DIR = "C:\\Users\\xinle\\Downloads\\three_layer_constant_fluxes_linear_hr144_Qu7.0e-04_Qb0.0e+00_f1.0e-04_Nh256_Nz128_strong_wind"

frame = 289
OUTPUT_PATH = "final_results\\LES_FC_$(frame).png"

# file = jldopen(joinpath(FILE_DIR, "xz_slice.jld2"))

# close(file)
plot_LES_3D(frame, FILE_DIR, OUTPUT_PATH, axis_images, title="", colorscheme=:turbo, rev=false)


text: include("Adam.jl")

max_stars_repo_name: JuliaNLSolvers/NLSolvers

max_issues_repo_name: JuliaNLSolvers/NLSolvers

max_forks_repo_name: JuliaNLSolvers/NLSolvers


text:
# Version management

"""
    version()

Returns the CUDA version as reported by the driver.
"""
@memoize function version()
    version_ref = Ref{Cint}()
    cuDriverGetVersion(version_ref)
    major, ver = divrem(version_ref[], 1000)
    minor, patch = divrem(ver, 10)
    return VersionNumber(major, minor, patch)
end

"""
    release()

Returns the CUDA release part of the version as returned by [`version`](@ref).
"""
@memoize release() = VersionNumber(version().major, version().minor)


text:
@testset "empty matrix" begin
    @test_throws MethodError BlockMatrix(0, 0)

    for A in [
        BlockMatrix(Array{AbstractArray{Int, 2}, 2}(undef, 0, 0)),
        BlockMatrix{Int}(Array{AbstractArray{Int, 2}, 2}(undef, 0, 0)),
        BlockMatrix{Int}(0, 0)
    ]
        @test A isa AbstractArray{Int, 2}
        @test size(A) == (0, 0)
        @test length(A) == 0
        @test_throws BoundsError A[]
        @test_throws BoundsError A[] = 42
        @test_throws BoundsError A[1]
        @test_throws BoundsError A[1] = 42
    end

    @test_throws AssertionError BlockMatrix(0, 0, zeros(Int, 2, 2))
    @test_throws AssertionError BlockMatrix{Int}(0, 0, zeros(Int, 2, 2))
    @test_throws MethodError BlockMatrix{Int}(0, 0, zeros(2, 2))
end

@testset "no rows" begin
    for A in [
        BlockMatrix(Array{AbstractArray{Int, 2}, 2}(undef, 0, 2)),
        BlockMatrix{Int}(Array{AbstractArray{Int, 2}, 2}(undef, 0, 2)),
        BlockMatrix{Int}(0, 2)
    ]
        @test A isa AbstractArray{Int, 2}
        @test size(A) == (0, 2)
        @test length(A) == 0
        @test_throws BoundsError A[]
        @test_throws BoundsError A[] = 42
        @test_throws BoundsError A[1]
        @test_throws BoundsError A[1] = 42
    end

    @test_throws AssertionError BlockMatrix(0, 2, zeros(Int, 2, 2))
    @test_throws AssertionError BlockMatrix{Int}(0, 2, zeros(Int, 2, 2))
    @test_throws MethodError BlockMatrix{Int}(0, 2, zeros(2, 2))
end

@testset "no cols" begin
    for A in [
        BlockMatrix(Array{AbstractArray{Int, 2}, 2}(undef, 2, 0)),
        BlockMatrix{Int}(Array{AbstractArray{Int, 2}, 2}(undef, 2, 0)),
        BlockMatrix{Int}(2, 0)
    ]
        @test A isa AbstractArray{Int, 2}
        @test size(A) == (2, 0)
        @test length(A) == 0
        @test_throws BoundsError A[]
        @test_throws BoundsError A[] = 42
        @test_throws BoundsError A[1]
        @test_throws BoundsError A[1] = 42
    end

    @test_throws AssertionError BlockMatrix(2, 0, zeros(Int, 2, 2))
    @test_throws AssertionError BlockMatrix{Int}(2, 0, zeros(Int, 2, 2))
    @test_throws MethodError BlockMatrix{Int}(2, 0, zeros(2, 2))
end

@testset "single element" begin
    B11 = 42 .* ones(Int, 1, 1)
    B = Array{AbstractArray{Int, 2}, 2}(undef, 1, 1)
    B[1] = B11

    for (i, A) in enumerate([
        BlockMatrix(B),
        BlockMatrix{Int}(B),
        BlockMatrix(1, 1, B11),
        BlockMatrix{Int}(1, 1, B11)
    ])
        @test A isa AbstractArray{Int, 2}
        @test size(A) == (1, 1)
        @test length(A) == 1
        A[] = 42 + i
        @test A[] isa Int
        @test A[] == B11[] == 42 + i
        A[1] = 43 + i
        @test A[1] isa Int
        @test A[1] == B11[1] == 43 + i
        A[1, 1] = 44 + i
        @test A[1, 1] isa Int
        @test A[1, 1] == B11[1, 1] == 44 + i
        @test_throws BoundsError A[0]
        @test_throws BoundsError A[0, 0]
        @test_throws BoundsError A[0, 1]
        @test_throws BoundsError A[1, 0]
        @test_throws BoundsError A[2]
        @test_throws BoundsError A[1, 2]
        @test_throws BoundsError A[2, 1]
        @test_throws BoundsError A[2, 2]
        @test_throws BoundsError A[0] = 42
        @test_throws BoundsError A[0, 0] = 42
        @test_throws BoundsError A[0, 1] = 42
        @test_throws BoundsError A[1, 0] = 42
        @test_throws BoundsError A[2] = 42
        @test_throws BoundsError A[1, 2] = 42
        @test_throws BoundsError A[2, 1] = 42
        @test_throws BoundsError A[2, 2] = 42
    end

    @test_throws AssertionError BlockMatrix{Int}(1, 1)
    @test_throws AssertionError BlockMatrix(1, 1, B11, B11)
    @test_throws AssertionError BlockMatrix{Int}(1, 1, B11, B11)
    @test_throws MethodError BlockMatrix{Float64}(1, 1, B11)
end

@testset "single row" begin
    B11 = [1 2]
    B12 = [3 4]
    B = Array{AbstractArray{Int, 2}, 2}(undef, 1, 2)
    B[1] = B11
    B[2] = B12

    for (i, A) in enumerate([
        BlockMatrix(B),
        BlockMatrix{Int}(B),
        BlockMatrix(1, 2, B11, B12),
        BlockMatrix{Int}(1, 2, B11, B12)
    ])
        @test A isa AbstractArray{Int, 2}
        @test size(A) == (1, 4)
        @test length(A) == 4
        for j in 1:4
            A[j] = 42i + j
            @test A[j] isa Int
            @test A[j] == 42i + j
            A[1, j] = 42i + j + 1
            @test A[1, j] isa Int
            @test A[1, j] == 42i + j + 1
        end
        @test_throws BoundsError A[]
        @test_throws BoundsError A[0]
        @test_throws BoundsError A[0, 0]
        @test_throws BoundsError A[0, 1]
        @test_throws BoundsError A[1, 0]
        @test_throws BoundsError A[5]
        @test_throws BoundsError A[1, 5]
        @test_throws BoundsError A[2, 1]
        @test_throws BoundsError A[] = 42
        @test_throws BoundsError A[0] = 42
        @test_throws BoundsError A[0, 0] = 42
        @test_throws BoundsError A[0, 1] = 42
        @test_throws BoundsError A[1, 0] = 42
        @test_throws BoundsError A[5] = 42
        @test_throws BoundsError A[1, 5] = 42
        @test_throws BoundsError A[2, 1] = 42
    end

    @test_throws AssertionError BlockMatrix{Int}(1, 2)
    @test_throws AssertionError BlockMatrix(1, 2, B11)
    @test_throws AssertionError BlockMatrix{Int}(1, 2, B11)
    @test_throws AssertionError BlockMatrix(1, 2, B11, B12, B12)
    @test_throws AssertionError BlockMatrix{Int}(1, 2, B11, B12, B12)
    @test_throws MethodError BlockMatrix{Float64}(1, 2, B11)
end

@testset "single col" begin
    B11 = reshape([1, 2], (2, 1))
    B21 = reshape([3, 4], (2, 1))
    B = Array{AbstractArray{Int, 2}, 2}(undef, 2, 1)
    B[1] = B11
    B[2] = B21

    for (i, A) in enumerate([
        BlockMatrix(B),
        BlockMatrix{Int}(B),
        BlockMatrix(2, 1, B11, B21),
        BlockMatrix{Int}(2, 1, B11, B21)
    ])
        @test A isa AbstractArray{Int, 2}
        @test size(A) == (4, 1)
        @test length(A) == 4
        for j in 1:4
            A[j] = 42i + j
            @test A[j] isa Int
            @test A[j] == 42i + j
            A[j, 1] = 42i + j + 1
            @test A[j, 1] isa Int
            @test A[j, 1] == 42i + j + 1
        end
        @test_throws BoundsError A[]
        @test_throws BoundsError A[0]
        @test_throws BoundsError A[0, 0]
        @test_throws BoundsError A[0, 1]
        @test_throws BoundsError A[1, 0]
        @test_throws BoundsError A[5]
        @test_throws BoundsError A[5, 1]
        @test_throws BoundsError A[1, 2]
        @test_throws BoundsError A[] = 42
        @test_throws BoundsError A[0] = 42
        @test_throws BoundsError A[0, 0] = 42
        @test_throws BoundsError A[0, 1] = 42
        @test_throws BoundsError A[1, 0] = 42
        @test_throws BoundsError A[5] = 42
        @test_throws BoundsError A[5, 1] = 42
        @test_throws BoundsError A[1, 2] = 42
    end

    @test_throws AssertionError BlockMatrix{Int}(2, 1)
    @test_throws AssertionError BlockMatrix(2, 1, B11)
    @test_throws AssertionError BlockMatrix{Int}(2, 1, B11)
    @test_throws AssertionError BlockMatrix(2, 1, B11, B21, B21)
    @test_throws AssertionError BlockMatrix{Int}(2, 1, B11, B21, B21)
    @test_throws MethodError BlockMatrix{Float64}(2, 1, B11)
end

@testset "single block" begin
    B11 = [1 2 3; 4 5 6]
    B = Array{AbstractArray{Int, 2}, 2}(undef, 1, 1)
    B[1] = B11

    for (i, A) in enumerate([
        BlockMatrix(B),
        BlockMatrix{Int}(B),
        BlockMatrix(1, 1, B11),
        BlockMatrix{Int}(1, 1, B11)
    ])
        @test A isa AbstractArray{Int, 2}
        @test size(A) == (2, 3)
        @test length(A) == 6
        for row in 1:2
            for col in 1:3
                A[row * col] = 42i + row * col
                @test A[row * col] isa Int
                @test A[row * col] == 42i + row * col
                A[row, col] = 42i + row * col + 1
                @test A[row, col] isa Int
                @test A[row, col] == 42i + row * col + 1
            end
        end
        @test_throws BoundsError A[]
        @test_throws BoundsError A[0]
        @test_throws BoundsError A[0, 0]
        @test_throws BoundsError A[0, 1]
        @test_throws BoundsError A[1, 0]
        @test_throws BoundsError A[7]
        @test_throws BoundsError A[1, 4]
        @test_throws BoundsError A[3, 1]
        @test_throws BoundsError A[] = 42
        @test_throws BoundsError A[0] = 42
        @test_throws BoundsError A[0, 0] = 42
        @test_throws BoundsError A[0, 1] = 42
        @test_throws BoundsError A[1, 0] = 42
        @test_throws BoundsError A[7] = 42
        @test_throws BoundsError A[1, 4] = 42
        @test_throws BoundsError A[3, 1] = 42
    end

    @test_throws AssertionError BlockMatrix{Int}(1, 1)
    @test_throws AssertionError BlockMatrix(1, 1, B11, B11)
    @test_throws AssertionError BlockMatrix{Int}(1, 1, B11, B11)
    @test_throws MethodError BlockMatrix{Float64}(1, 1, B11)
end

@testset "single block row" begin
    B11 = [1 2 3; 7 8 9]
    B12 = [4 5 6; 10 11 12]
    B = Array{AbstractArray{Int, 2}, 2}(undef, 1, 2)
    B[1] = B11
    B[2] = B12

    for (i, A) in enumerate([
        BlockMatrix(B),
        BlockMatrix{Int}(B),
        BlockMatrix(1, 2, B11, B12),
        BlockMatrix{Int}(1, 2, B11, B12)
    ])
        @test A isa AbstractArray{Int, 2}
        @test size(A) == (2, 6)
        @test length(A) == 12
        for row in 1:2
            for col in 1:6
                A[row * col] = 42i + row * col
                @test A[row * col] isa Int
                @test A[row * col] == 42i + row * col
                A[row, col] = 42i + row * col + 1
                @test A[row, col] isa Int
                @test A[row, col] == 42i + row * col + 1
            end
        end
        @test_throws BoundsError A[]
        @test_throws BoundsError A[0]
        @test_throws BoundsError A[0, 0]
        @test_throws BoundsError A[0, 1]
        @test_throws BoundsError A[1, 0]
        @test_throws BoundsError A[13]
        @test_throws BoundsError A[1, 7]
        @test_throws BoundsError A[3, 1]
        @test_throws BoundsError A[] = 42
        @test_throws BoundsError A[0] = 42
        @test_throws BoundsError A[0, 0] = 42
        @test_throws BoundsError A[0, 1] = 42
        @test_throws BoundsError A[1, 0] = 42
        @test_throws BoundsError A[13] = 42
        @test_throws BoundsError A[1, 7] = 42
        @test_throws BoundsError A[3, 1] = 42
    end

    @test_throws AssertionError BlockMatrix(1, 2, B11)
    @test_throws AssertionError BlockMatrix{Int}(1, 2, B11)
    @test_throws AssertionError BlockMatrix(1, 2, B11, B12, B12)
    @test_throws AssertionError BlockMatrix{Int}(1, 2, B11, B12, B12)
    @test_throws MethodError BlockMatrix{Float64}(1, 2, B11, B12)
end

@testset "single block col" begin
    B11 = [1 2 3; 4 5 6]
    B21 = [7 8 9; 10 11 12]
    B = Array{AbstractArray{Int, 2}, 2}(undef, 2, 1)
    B[1] = B11
    B[2] = B21

    for (i, A) in enumerate([
        BlockMatrix(B),
        BlockMatrix{Int}(B),
        BlockMatrix(2, 1, B11, B21),
        BlockMatrix{Int}(2, 1, B11, B21)
    ])
        @test A isa AbstractArray{Int, 2}
        @test size(A) == (4, 3)
        @test length(A) == 12
        for row in 1:4
            for col in 1:3
                A[row * col] = 42i + row * col
                @test A[row * col] isa Int
                @test A[row * col] == 42i + row * col
                A[row, col] = 42i + row * col + 1
                @test A[row, col] isa Int
                @test A[row, col] == 42i + row * col + 1
            end
        end
        @test_throws BoundsError A[]
        @test_throws BoundsError A[0]
        @test_throws BoundsError A[0, 0]
        @test_throws BoundsError A[0, 1]
        @test_throws BoundsError A[1, 0]
        @test_throws BoundsError A[13]
        @test_throws BoundsError A[1, 4]
        @test_throws BoundsError A[5, 1]
        @test_throws BoundsError A[] = 42
        @test_throws BoundsError A[0] = 42
        @test_throws BoundsError A[0, 0] = 42
        @test_throws BoundsError A[0, 1] = 42
        @test_throws BoundsError A[1, 0] = 42
        @test_throws BoundsError A[13] = 42
        @test_throws BoundsError A[1, 4] = 42
        @test_throws BoundsError A[5, 1] = 42
    end

    @test_throws AssertionError BlockMatrix(2, 1, B11)
    @test_throws AssertionError BlockMatrix{Int}(2, 1, B11)
    @test_throws AssertionError BlockMatrix(2, 1, B11, B21, B21)
    @test_throws AssertionError BlockMatrix{Int}(2, 1, B11, B21, B21)
    @test_throws MethodError BlockMatrix{Float64}(2, 1, B11, B21)
end

@testset "ordinary" begin
    B11 = [1 2 3; 7 8 9]
    B12 = [4 5 6; 10 11 12]
    B21 = [13 14 15; 19 20 21]
    B22 = [16 17 18; 22 23 24]
    B = Array{AbstractArray{Int, 2}, 2}(undef, 2, 2)
    B[1] = B11
    B[2] = B21
    B[3] = B12
    B[4] = B22

    for (i, A) in enumerate([
        BlockMatrix(B),
        BlockMatrix{Int}(B),
        BlockMatrix(2, 2, B11, B12, B21, B22),
        BlockMatrix{Int}(2, 2, B11, B12, B21, B22)
    ])
        @test A isa AbstractArray{Int, 2}
        @test size(A) == (4, 6)
        @test length(A) == 24
        for row in 1:4
            for col in 1:6
                A[row * col] = 42i + row * col
                @test A[row * col] isa Int
                @test A[row * col] == 42i + row * col
                A[row, col] = 42i + row * col + 1
                @test A[row, col] isa Int
                @test A[row, col] == 42i + row * col + 1
            end
        end
        @test_throws BoundsError A[]
        @test_throws BoundsError A[0]
        @test_throws BoundsError A[0, 0]
        @test_throws BoundsError A[0, 1]
        @test_throws BoundsError A[1, 0]
        @test_throws BoundsError A[25]
        @test_throws BoundsError A[1, 7]
        @test_throws BoundsError A[5, 1]
        @test_throws BoundsError A[] = 42
        @test_throws BoundsError A[0] = 42
        @test_throws BoundsError A[0, 0] = 42
        @test_throws BoundsError A[0, 1] = 42
        @test_throws BoundsError A[1, 0] = 42
        @test_throws BoundsError A[25] = 42
        @test_throws BoundsError A[1, 7] = 42
        @test_throws BoundsError A[5, 1] = 42
    end

    @test_throws AssertionError BlockMatrix(2, 2, B11, B12, B21)
    @test_throws AssertionError BlockMatrix{Int}(2, 2, B11, B12, B21)
    @test_throws AssertionError BlockMatrix(2, 2, B11, B12, B21, B22, B22)
    @test_throws AssertionError BlockMatrix{Int}(2, 2, B11, B12, B21, B22, B22)
    @test_throws MethodError BlockMatrix{Float64}(2, 2, B11, B12, B21, B22)
end


text:
export gauss_legendre_weights;   # For debugging only
export distributed_kernel_trapezoidal;    # For debugging only
export distributed_kernel_gauss_legendre; # For debugging only

"
#Creates the NEP associated with example in E. Jarlebring and
#W. Michiels and K. Meerbergen, The infinite  {Arnoldi} method and an application to time-delay systems with distributed delays, Delay Systems - Methods, Applications and New Trends, 2012
"
function gallery_dep_distributed()
    # E. Jarlebring and W. Michiels and K. Meerbergen},
    # The infinite  {Arnoldi} method and an application to time-delay systems with distributed delays,
    # Delay Systems - Methods, Applications and New Trends
    # 2012
    #
    # Some correct eigenvalues:
    # -0.400236388049641 + 0.970633098237807i
    # -0.400236388049641 - 0.970633098237807i
    #  2.726146249832675 + 0.000000000000000i
    # -1.955643591177653 + 3.364550574688863i
    # -1.955643591177653 - 3.364550574688863i
    #  4.493937056300693 + 0.000000000000000i
    # -1.631513006819252 + 4.555484848248613i
    # -1.631513006819252 - 4.555484848248613i
    # -1.677320660400946 + 7.496870451838560i
    # -1.677320660400946 - 7.496870451838560i
    #
    #

    A0 = Matrix(-1.0I, 3, 3)
    A1=[2.5    2.8   -0.5
        1.8    0.3    0.3
        -2.3   -1.4    3.5];
    A2=[1.7    0.7   -0.3
        -2.4   -2.1   -0.2
        2.0    0.7    0.4];
    A3=[1.4   -1.3    0.4
        1.4    0.7    1.0
        0.6    1.6    1.7];
    idop= S -> S
    oneop = S -> one(S)
    f1 = S -> exp(-S)
    N = 10
    f2 = S -> distributed_kernel_gauss_legendre(S, N)
    #N = 1000
    #f2 = S -> distributed_kernel_trapezoidal(Matrix(S), N)
    return SPMF_NEP([A0,A1,A2,A3],[idop,oneop,f1,f2])
end

function distributed_kernel_trapezoidal(S,N0)

    fS = x -> exp(x*S) * (exp((x+0.5)^2) - exp(1/4))

    F=zero(S)

    h=1/N0;
    for i=1:(N0+1) # Trapezoidal rule
        x=float((i-1)*h-1)
        if (i==1 || i==N0+1)
            F += 0.5*fS(x)*h;
        else
            F += fS(x)*h;
        end
    end
    return F
end
"""
Computes distributed kernel matrix function using
Gauss-Legender quadrature.
"""
function distributed_kernel_gauss_legendre(S,N)
    f=x-> (exp((x+0.5)^2)-exp(1/4))
    F=zero(S);
    xv,wv=gauss_legendre_weights(N,-1,0);
    local E = one(S)

    accumulative_expm_comp=true
    for i=1:length(xv)
        if (accumulative_expm_comp)
            # An accumulative way to compute E=exp(xv[i]*S) which is
            # faster due to the fact that scaling and squaring
            # for exp((xv[i]-xv[i-1])*S) is faster than
            # exp(xv[i]*S)
            if (i==1)
                E=exp(xv[1]*S);
            else
                E=E*exp((xv[i]-xv[i-1])*S)
            end
        else
            E=exp(xv[i]*S);

        end

        fSw=E*(f(xv[i])*wv[i]);
        F=F+fSw;
    end
    return F
end

function  gauss_legendre_weights(N,a,b)

    N1=N; N2=N+1;
    xu = range(-1, stop = 1, length = N1)

    # L will be the Legendre-Gauss Vandermonde Matrix
    L=zeros(N1,N2);

    # "Derivative" of L
    Lp=zeros(N1,N2);

    # Starting values of Newton's method
    y = cos.((2 * (0:(N-1)) .+ 1) * pi / (2 * (N-1) + 2)) + (0.27/N1) * sin.(pi * xu * (N-1) / N2)
    y0 = 2

    local Lp0
    # Newton's method to decide points
    # Iterate until new points are uniformly within epsilon of old points

    while maximum(abs.(y .- y0)) > eps()

        L[:,1]=ones(size(L,1));
        Lp[:,1]=zeros(size(Lp,1));

        L[:,2]=y;
        Lp[:,2]=ones(size(Lp,1));

        # Recurrence
        for k=2:N1
            L[:,k+1]=( (2*k-1)*(y).*L[:,k]-(k-1)*L[:,k-1] )/k;
        end

        Lp0 = N2 * (L[:,N1] - y .* L[:,N2]) ./ (1 .- y.^2)

        y0=y;
        y=y0-L[:,N2]./Lp0;
    end

    # Linear map from[-1,1] to [a,b]
    x = (a * (1 .- y) + b * (1 .+ y)) / 2

    # Compute the weights
    w = (b-a) ./ ((1 .- y.^2) .* Lp0.^2) * (N2/N1)^2

    return x,w
end


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

include("ahmc_integrator.jl")
include("ahmc_proposal.jl")
include("ahmc_metric.jl")
include("ahmc_sampleid.jl")
include("ahmc_sampler.jl")
include("ahmc_adaptor.jl")


text:
using MinimallyDisruptiveCurves
using Test

@testset "mass_spring" begin
            include("mass_spring.jl")
        end


text:
function create(c::Client, ::Type{Webhook}, ch::DiscordChannel; kwargs...)
    return create_webhook(c, ch.id; kwargs...)
end

retrieve(c::Client, ::Type{Webhook}, ch::DiscordChannel) = get_channel_webhooks(c, ch.id)
retrieve(c::Client, ::Type{Webhook}, g::AbstractGuild) = get_guild_webhooks(c, g.id)
retrieve(c::Client, ::Type{Webhook}, webhook::Integer) = get_webhook(c, webhook)
function retrieve(c::Client, ::Type{Webhook}, webhook::Integer, token::AbstractString)
    return get_webhook_with_token(c, webhook, token)
end

update(c::Client, w::Webhook; kwargs...) = modify_webhook(c, w.id; kwargs...)
function update(c::Client, w::Webhook, token::AbstractString; kwargs...)
    return modify_webhook_with_token(c, w.id, token; kwargs...)
end

delete(c::Client, w::Webhook) = delete_webhook(c, w.id)
function delete(c::Client, w::Webhook, token::AbstractString)
    return delete_webhook_with_token(c, w.id, token)
end


text:
@reexport using DataFrames

export sourceParticles

foreach([<=, >=, <, >]) do cmp
    MonteCarloMeasurements.register_primitive(cmp, eval)
end

export makeParticles
function makeParticles(m :: Model)
    fpre = @eval $(sourceParticles(m))
    f(;kwargs...) = Base.invokelatest(fpre; kwargs...)
end

export particles
particles(m::Model; kwargs...) = makeParticles(m)(;kwargs...)

function particles(m::Model, n::Int64; kwargs...)
    r = makeParticles(m)
    [r(;kwargs...) for j in 1:n] |> DataFrame

end

particles(d) = Particles(1000, d)

function sourceParticles(m::Model)
    m = canonical(m)
    proc(m, st::Let)     = :($(st.x) = $(st.rhs))
    proc(m, st::Follows) = :($(st.x) = parts($(st.rhs)))
    proc(m, st::Return)  = :(return $(st.rhs))
    proc(m, st::LineNumber) = nothing

    body = buildSource(m, proc) |> striplines

    argsExpr = Expr(:tuple,freeVariables(m)...)

    stochExpr = begin
        vals = map(variables(m)) do x Expr(:(=), x,x) end
        Expr(:tuple, vals...)
    end

    @gensym particles

    flatten(@q (
        function $particles(args...;kwargs...)
            @unpack $argsExpr = kwargs
            $body
            $stochExpr
        end
    ))

end

particles(n :: NUTS_result) = particles(n.samples)

function particles(s::Vector{NamedTuple{vs, T}} where {vs, T})
    nt = NamedTuple()
    for v in keys(s[1])
        nt = merge(nt, [v => Particles(getproperty.(s,v))])
    end
    nt
end

export parts

# Just a little helper function for particles
# https://github.com/baggepinnen/MonteCarloMeasurements.jl/issues/22
parts(m::Model; N=1000) = particles(m)
parts(x::Normal{P} where {P <: AbstractParticles}; N=1000) = Particles(length(x.μ), Normal()) * x.σ + x.μ
parts(x::Sampleable{F,S}; N=1000) where {F,S} = Particles(N,x)
parts(x::Integer; N=1000) = parts(float(x))
parts(x::Real; N=1000) = parts(repeat([x],N))
parts(x::AbstractArray; N=1000) = Particles(x)
parts(p::Particles; N=1000) = p
parts(d::For; N=1000) = map(d.θs) do θ
    parts(d.f(θ))
end

parts(d::iid; N=1000) = map(1:d.size) do j parts(d.dist) end
# size
# dist

# MonteCarloMeasurements.Particles(n::Int, d::For) =
# -(a::Particles{Float64,1000}, b::Array{Float64,1}) = [a-bj for bj in b]

# Base.promote(a::Particles{T,N}, b)  where {T,N} = (a,parts(b))
# Base.promote(a, b::Particles{T,N})  where {T,N} = (parts(a),b)

# promote_rule(::Type{A}, ::Type{B}) where {A <: Real, B <: AbstractParticles{T,N}} where {T} = AbsractParticles{promote_type(A,T),N} where {N}
# promote_rule(::Type{B}, ::Type{A}) where {A <: Real, B <: AbstractParticles{T,N}} where {T} = AbsractParticles{promote_type(A,T),N} where {N}


text:
include("StructEquality.jl")

a_counter = CounterX(1)

max_stars_repo_name: Sohva/oolong

max_issues_repo_name: Sohva/oolong

max_forks_repo_name: oolongtensor/oolong


text:
using SolarLEDs
using Test


text:
julia> decode(Float32(1/3))
('0', "01111111", "10101010101010101010101")

julia> decode(Float64(1/3))
('0', "01111111", "1010101010101010101010101010101010101010101010101010101")

julia> decode(Float64(Float32(1/3)))
('0', "01111111", "1010101010101010101010101100000000000000000000000000000")


text:
using GraphTools
using Base.Test

u = @graph
@test "graph graphname {\n}\n" == to_dot(u)

u = @graph A
@test "graph graphname {\n1 [\"label\"=\"A\"]\n}\n" == to_dot(u)

u = @graph A - B
@test "graph graphname {\n1 [\"label\"=\"A\"]\n1 -- 2 \n2 [\"label\"=\"B\"]\n}\n" == to_dot(u)

u = @graph A - B - C
@test "graph graphname {\n1 [\"label\"=\"A\"]\n1 -- 2 \n2 [\"label\"=\"B\"]\n2 -- 3 \n3 [\"label\"=\"C\"]\n}\n" == to_dot(u)

u = @graph A - B - C - A
@test "graph graphname {\n1 [\"label\"=\"A\"]\n1 -- 2 \n1 -- 3 \n2 [\"label\"=\"B\"]\n2 -- 3 \n3 [\"label\"=\"C\"]\n}\n" == to_dot(u)
@test 1 == indexof(u, "A")
@test 2 == indexof(u, "B")
@test 0 == indexof(u, "Z")

u = @graph A 50- B 60- C
@test "graph graphname {\n1 [\"label\"=\"A\"]\n1 -- 2 [\"label\"=\"50\",\"weight\"=\"50\"]\n2 [\"label\"=\"B\"]\n2 -- 3 [\"label\"=\"60\",\"weight\"=\"60\"]\n3 [\"label\"=\"C\"]\n}\n" == to_dot(u)

u = @graph A 0.5- B 0.6- C
@test "graph graphname {\n1 [\"label\"=\"A\"]\n1 -- 2 [\"label\"=\"0.5\",\"weight\"=\"0.5\"]\n2 [\"label\"=\"B\"]\n2 -- 3 [\"label\"=\"0.6\",\"weight\"=\"0.6\"]\n3 [\"label\"=\"C\"]\n}\n" == to_dot(u)

ua = @graph A
ub = @graph B
@test "graph graphname {\n1 [\"label\"=\"A\"]\n2 [\"label\"=\"B\"]\n}\n" == to_dot(ua + ub)

u = @graph "A"
@test "graph graphname {\n1 [\"label\"=\"A\"]\n}\n" == to_dot(u)

u = @graph "A" - "B"
@test "graph graphname {\n1 [\"label\"=\"A\"]\n1 -- 2 \n2 [\"label\"=\"B\"]\n}\n" == to_dot(u)

u = @graph "A" - "B" - "C"
@test "graph graphname {\n1 [\"label\"=\"A\"]\n1 -- 2 \n2 [\"label\"=\"B\"]\n2 -- 3 \n3 [\"label\"=\"C\"]\n}\n" == to_dot(u)

u = @graph "A" - "B" - "C" - "A"
@test "graph graphname {\n1 [\"label\"=\"A\"]\n1 -- 2 \n1 -- 3 \n2 [\"label\"=\"B\"]\n2 -- 3 \n3 [\"label\"=\"C\"]\n}\n" == to_dot(u)
@test 1 == indexof(u, "A")
@test 2 == indexof(u, "B")
@test 0 == indexof(u, "Z")

u = @graph "A" 50- "B" 60- "C"
@test "graph graphname {\n1 [\"label\"=\"A\"]\n1 -- 2 [\"label\"=\"50\",\"weight\"=\"50\"]\n2 [\"label\"=\"B\"]\n2 -- 3 [\"label\"=\"60\",\"weight\"=\"60\"]\n3 [\"label\"=\"C\"]\n}\n" == to_dot(u)

u = @graph "A" 0.5- "B" 0.6- "C"
@test "graph graphname {\n1 [\"label\"=\"A\"]\n1 -- 2 [\"label\"=\"0.5\",\"weight\"=\"0.5\"]\n2 [\"label\"=\"B\"]\n2 -- 3 [\"label\"=\"0.6\",\"weight\"=\"0.6\"]\n3 [\"label\"=\"C\"]\n}\n" == to_dot(u)

ua = @graph "A"
ub = @graph "B"
@test "graph graphname {\n1 [\"label\"=\"A\"]\n2 [\"label\"=\"B\"]\n}\n" == to_dot(ua + ub)


text:
using TensorKit
using TensorOperations
using LinearAlgebra
using statmech_tm_solver

io = open("result_triangular_ising_fulldiag_alternative.txt", "w")
close(io)
T0 = mpo_triangular_AF_ising_alternative()
T1, T2 = mpo_triangular_AF_ising_adapter()
t_fuse = isomorphism(ℂ^32, ℂ^4*ℂ^4*ℂ^2)

@tensor T[-1, -2; -3, -4] := t_fuse[-1, 3, 4, 5] * T1[3, -2, 1, 6] * T0[4, 1, 2, 7] * T2[5, 2, -3, 8] * t_fuse'[6, 7, 8, -4]

U, S, V, ϵ = tsvd(permute(T, (1, 2, 3), (4,)), trunc=truncerr(1e-10))
T = permute(S*V*permute(U, (1, ), (2, 3, 4)), (1, 2), (3, 4))

T5 = ncon([T, T, T, T, T], [[1, -1, -2, 2], [2, -3, -4, 3], [3, -5, -6, 4], [4, -7, -8, 5], [5, -9, -10, 1]])
T5 = permute(T5, (1, 3, 5, 7, 9), (2, 4, 6, 8, 10))

T6 = ncon([T, T, T, T, T, T], [[1, -1, -2, 2], [2, -3, -4, 3], [3, -5, -6, 4], [4, -7, -8, 5], [5, -9, -10, 6], [6, -11, -12, 1]])
T6 = permute(T6, (1, 3, 5, 7, 9, 11), (2, 4, 6, 8, 10, 12))

T7 = ncon([T, T, T, T, T, T, T], [[1, -1, -2, 2], [2, -3, -4, 3], [3, -5, -6, 4], [4, -7, -8, 5], [5, -9, -10, 6], [6, -11, -12, 7], [7, -13, -14, 1]])
T7 = permute(T7, (1, 3, 5, 7, 9, 11, 13), (2, 4, 6, 8, 10, 12, 14))

T8 = ncon([T, T, T, T, T, T, T, T], [[1, -1, -2, 2], [2, -3, -4, 3], [3, -5, -6, 4], [4, -7, -8, 5], [5, -9, -10, 6], [6, -11, -12, 7], [7, -13, -14, 8], [8, -15, -16, 1]])
T8 = permute(T8, (1, 3, 5, 7, 9, 11, 13, 15), (2, 4, 6, 8, 10, 12, 14, 16))

T9 = ncon([T, T, T, T, T, T, T, T, T], [[1, -1, -2, 2], [2, -3, -4, 3], [3, -5, -6, 4], [4, -7, -8, 5], [5, -9, -10, 6], [6, -11, -12, 7], [7, -13, -14, 8], [8, -15, -16, 9], [9, -17, -18, 1]])
T9 = permute(T9, (1, 3, 5, 7, 9, 11, 13, 15, 17), (2, 4, 6, 8, 10, 12, 14, 16, 18))

T10 = ncon([T, T, T, T, T, T, T, T, T, T], [[1, -1, -2, 2], [2, -3, -4, 3], [3, -5, -6, 4], [4, -7, -8, 5], [5, -9, -10, 6], [6, -11, -12, 7], [7, -13, -14, 8], [8, -15, -16, 9], [9, -17, -18, 10], [10, -19, -20, 1]])
T10 = permute(T10, (1, 3, 5, 7, 9, 11, 13, 15, 17, 19), (2, 4, 6, 8, 10, 12, 14, 16, 18, 20))

T11 = ncon([T, T, T, T, T, T, T, T, T, T, T], [[1, -1, -2, 2], [2, -3, -4, 3], [3, -5, -6, 4], [4, -7, -8, 5], [5, -9, -10, 6], [6, -11, -12, 7], [7, -13, -14, 8], [8, -15, -16, 9], [9, -17, -18, 10], [10, -19, -20, 11], [11, -21, -22, 1]])
T11 = permute(T11, (1, 3, 5, 7, 9, 11, 13, 15, 17, 19, 21), (2, 4, 6, 8, 10, 12, 14, 16, 18, 20, 22))

T12 = ncon([T, T, T, T, T, T, T, T, T, T, T, T], [[1, -1, -2, 2], [2, -3, -4, 3], [3, -5, -6, 4], [4, -7, -8, 5], [5, -9, -10, 6], [6, -11, -12, 7], [7, -13, -14, 8], [8, -15, -16, 9], [9, -17, -18, 10], [10, -19, -20, 11], [11, -21, -22, 12], [12, -23, -24, 1]])
T12 = permute(T12, (1, 3, 5, 7, 9, 11, 13, 15, 17, 19, 21, 23), (2, 4, 6, 8, 10, 12, 14, 16, 18, 20, 22, 24))

for (n, Tn) in zip([5, 6, 7, 8, 9, 10, 11, 12], [T5, T6, T7, T8, T9, T10, T11, T12])
    w, _ = eig(Tn)
    w = diag(w.data)
    w_reals = real.(w)
    w_imags = imag.(w)

    io = open("result_triangular_ising_fulldiag_alternative.txt", "a+")
    for (w_real, w_imag) in zip(w_reals, w_imags)
        println(io, "$(n)   $(w_real)   $(w_imag) ")
    end
    close(io)
end


text:
# ---
# title: 83. Remove Duplicates from Sorted List
# id: problem83
# author: Tian Jun
# date: 2020-10-31
# difficulty: Easy
# categories: Linked List
# link: <https://leetcode.com/problems/remove-duplicates-from-sorted-list/description/>
# hidden: true
# ---
#
# Given a sorted linked list, delete all duplicates such that each element
# appear only _once_.
#
# **Example 1:**
#
#
#
#     Input: 1->1->2
#     Output: 1->2
#
#
# **Example 2:**
#
#
#
#     Input: 1->1->2->3->3
#     Output: 1->2->3
#
#
#
## @lc code=start
using LeetCode

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


text:
# # [Basic Epidemiology Models](@id epidemiology_example)
#
#md # [![](https://img.shields.io/badge/show-nbviewer-579ACA.svg)](@__NBVIEWER_ROOT_URL__/examples/epidemiology.ipynb)

using Petri
using LabelledArrays
using Plots
using DiffEqJump
using StochasticDiffEq
using OrdinaryDiffEq

# ### SIR Model
#
# The SIR model represents the epidemiological dynamics of an infectious disease
# that causes immunity in its victims. There are three *states:* `Suceptible
# ,Infected, Recovered`. These states interact through two *transitions*.
# Infection has the form `S+I -> 2I` where a susceptible person meets an
# infected person and results in two infected people. The second transition is
# recovery `I -> R` where an infected person recovers spontaneously.

S  = [:S,:I,:R]
Δ  = LVector(
       inf=(LVector(S=1, I=1), LVector(I=2)),
       rec=(LVector(I=1),      LVector(R=1)),
     )
sir = Petri.Model(S, Δ)

Graph(sir)

# Once a model is defined, we can define out initial parameters `u0`, a time
# span `tspan`, and the transition rates of the interactions `β`

u0 = LVector(S=990.0, I=10.0, R=0.0)
tspan = (0.0,40.0)
# add a dynamic transition rate for infection
# where the rate of infection decreases over time
# as is dependent on the current state of the system
β = LVector(inf=((u,t)->((3/sum(u))/(t+1))), rec=0.25);

# each transition rates can one of three options:
#
# - constant: `β = [.25]`
#   - where the rate is specified by a value of type `Number`
# - time dependent: `β = [t->((3/1000)/(t+1))]`
#   - where `t` is the current time step
# - state and time dependent: `β = [(u,t)->((3/sum(u))/(t+1))]`
#   - where `u` is the current state of `u0` and `t` is the current time step

# Petri.jl provides interfaces to StochasticDiffEq.jl, DiffEqJump.jl, and
# OrdinaryDiffEq.jl Here, we call the `JumpProblem` function that returns an
# DiffEqJump problem object that can be passed to the DiffEqJump solver which
# can then be plotted and visualized

prob = JumpProblem(sir, u0, tspan, β)
sol = DiffEqJump.solve(prob,SSAStepper())

plot(sol)

# Similarly, we can generated `SDEProblem` statements that can be used with
# StochasticDiffEq solvers

prob, cb = SDEProblem(sir, u0, tspan, β)
sol = StochasticDiffEq.solve(prob,LambaEM(),callback=cb)

plot(sol)

# Lastly, we can generated `ODEProblem` statements that can be used with
# OrdinOrdinaryDiffEq solvers

prob = ODEProblem(sir, u0, tspan, β)
sol = OrdinaryDiffEq.solve(prob,Tsit5(),reltol=1e-8,abstol=1e-8)

plot(sol)

# ### SEIR Model

S = [:S,:E,:I,:R]
Δ = LVector(
      exp=(LVector(S=1, I=1), LVector(I=1, E=1)),
      inf=(LVector(E=1),      LVector(I=1)),
      rec=(LVector(I=1),      LVector(R=1)),
    )
seir = Petri.Model(S, Δ)

Graph(seir)
#-
u0 = LVector(S=990.0, E=10.0, I=0.0, R=0.0)
tspan = (0.0,40.0)
β = LVector(exp=0.7/sum(u0), inf=0.5, rec=0.25)

prob, cb = SDEProblem(seir, u0, tspan, β)
sol = StochasticDiffEq.solve(prob,LambaEM(),callback=cb)

plot(sol)

# ### SEIRD Model

S = [:S,:E,:I,:R, :D]
Δ = LVector(
      exp=(LVector(S=1, I=1), LVector(I=1, E=1)),
      inf=(LVector(E=1),      LVector(I=1)),
      rec=(LVector(I=1),      LVector(R=1)),
      die=(LVector(I=1),      LVector(D=1)),
    )
seird = Petri.Model(S, Δ)

Graph(seird)
#-
u0 = LVector(S=990.0, E=10.0, I=0.0, R=0.0, D=0.0)
tspan = (0.0,40.0)
β = LVector(exp=0.9/sum(u0), inf=0.9, rec=0.25, die=0.03)

prob, cb = SDEProblem(seird, u0, tspan, β)
sol = StochasticDiffEq.solve(prob,LambaEM(),callback=cb)

plot(sol)


text:
## Example script for constant-strength source panel method
using LinearAlgebra
using Base.Iterators
using Seaborn
using AeroMDAO

## Airfoil
airfoil = Foil(naca4((0,0,1,2), 81; sharp_trailing_edge = true))
V, α    = 1., 0.
ρ       = 1.225
uniform = Uniform2D(V, α)
num_pans = 80
panels  = paneller(airfoil, num_pans);

## Constant-strength source panel
A = @time constant_source_matrix(panels)
b = @time constant_source_boundary_condition(panels, velocity(uniform))
σs = A \ b

##
sum(σs)

## Pressure coefficient
pts = collocation_point.(panels);
panel_vels = [ velocity(uniform) .+ sum(source_velocity.(σs, panels, Ref(panel_i))) for panel_i in panels ]

qts = @. dot(panel_vels, panel_tangent(panels))
cps = @. 1 - qts^2 / uniform.mag^2

##
upper, lower = get_surface_values(panels, cps)
lower = [ upper[end]; lower; upper[1] ]

plot(first.(upper), last.(upper), label = "Upper")
plot(first.(lower), last.(lower), label = "Lower")
ylim(maximum(cps), minimum(cps))
xlabel("(x/c)")
ylabel("Cp")
# legend()
show()

## Plotting
x_domain, y_domain = (-1, 2), (-1, 1)
grid_size = 50
x_dom, y_dom = range(x_domain..., length = grid_size), range(y_domain..., length = grid_size)
grid = product(x_dom, y_dom)
pts = panel_points(panels);

source_vels = [ velocity(uniform) .+ sum(source_velocity.(σs, panels, x, y)) for (x, y) in grid ]

speeds = @. sqrt(first(source_vels)^2 + last(source_vels)^2)
cps = @. 1 - speeds^2 / uniform.mag^2

contourf(first.(grid), last.(grid), cps)
CB1 = colorbar(label = "Pressure Coefficient (Cp)")
# quiver(first.(grid), last.(grid), first.(source_vels), last.(source_vels), speeds)
streamplot(first.(grid)', last.(grid)', first.(source_vels)', last.(source_vels)', color = speeds', cmap = "coolwarm", density = 3)
CB2 = colorbar(orientation="horizontal", label = "Relative Airspeed (U/U∞)")
fill(first.(pts), last.(pts), color = "black", zorder = 3)
tight_layout()
show()

max_forks_repo_name: GodotMisogi/AeroMDAO


text:
## Utilities and solutions for path search algorithms ##

mutable struct PathNode
    id::UInt
    state::State
    path_cost::Float64
    parent_id::Union{UInt,Nothing}
    parent_action::Union{Term,Nothing}
end

PathNode(id, state, path_cost) =
    PathNode(id, state, path_cost, nothing, nothing)

function reconstruct(node_id::UInt, search_tree::Dict{UInt,PathNode})
    plan, traj = Term[], State[search_tree[node_id].state]
    while node_id in keys(search_tree)
        node = search_tree[node_id]
        if node.parent_id === nothing break end
        pushfirst!(plan, node.parent_action)
        pushfirst!(traj, node.state)
        node_id = node.parent_id
    end
    return plan, traj
end

"Solution type for search-based planners that produce fully ordered plans."
mutable struct PathSearchSolution{T} <: OrderedSolution
    status::Symbol
    plan::Vector{Term}
    trajectory::Union{Vector{State},Nothing}
    search_tree::Union{Dict{UInt,PathNode},Nothing}
    frontier::T
end

PathSearchSolution(status::Symbol, plan) =
    PathSearchSolution(status, plan, nothing, nothing, nothing)
PathSearchSolution(status::Symbol, plan, trajectory) =
    PathSearchSolution(status, plan, trajectory, nothing, nothing)

get_action(sol::OrderedPlan, t::Int, state::State) = sol.plan[t]

Base.iterate(sol::PathSearchSolution) = iterate(sol.plan)
Base.iterate(sol::PathSearchSolution, istate) = iterate(sol.plan, istate)
Base.getindex(sol::PathSearchSolution, i::Int) = getindex(sol.plan, i)
Base.length(sol::PathSearchSolution) = length(sol.plan)


text:
#!/bin/env julia

const filename = "julia-vim-L2U-table"

# Keeping this from an old revision, just in case:
# Combining chars ranges, as obtained from https://en.wikipedia.org/wiki/Combining_character:
# ('\u0300' ≤ u ≤ '\u036F') ||
# ('\u1AB0' ≤ u ≤ '\u1AFF') ||
# ('\u1DC0' ≤ u ≤ '\u1DFF') ||
# ('\u20D0' ≤ u ≤ '\u20FF') ||
# ('\uFE20' ≤ u ≤ '\uFE2F')

# Most of this code is copy-pasted and slightly adapted from here:
# JULIA_HOME/doc/src/manual/unicode-input.md b/doc/src/manual/unicode-input.md
# except that instead of producing Markdown we want to generate a vim
# documentation file.

function tab_completions(symbols...)
    completions = Dict{String, Vector{String}}()
    for each in symbols, (k, v) in each
        completions[v] = push!(get!(completions, v, String[]), k)
    end
    return completions
end

function unicode_data()
    file = normpath(JULIA_HOME, "..", "..", "doc", "UnicodeData.txt")
    names = Dict{UInt32, String}()
    open(file) do unidata
        for line in readlines(unidata)
            id, name, desc = split(line, ";")[[1, 2, 11]]
            codepoint = parse(UInt32, "0x$id")
            names[codepoint] = (name == "" ? desc : desc == "" ? name : "$name / $desc")
        end
    end
    return names
end

# Prepend a dotted circle ('◌' i.e. '\u25CC') to combining characters
function fix_combining_chars(char)
    cat = Base.UTF8proc.category_code(char)
    return string(cat == 6 || cat == 8 ? "◌" : "", char)
end

function table_entries(completions, unicode_dict)
    code = String[]
    unicode = String[]
    latex = String[]
    desc = String[]

    for (chars, inputs) in sort!(collect(completions), by = first)
        code_points, unicode_names, characters = String[], String[], String[]
        for char in chars
            push!(code_points, "U+$(uppercase(hex(char, 5)))")
            push!(unicode_names, get(unicode_dict, UInt32(char), "(No Unicode name)"))
            push!(characters, isempty(characters) ? fix_combining_chars(char) : "$char")
        end
        push!(code, join(code_points, " + "))
        push!(unicode, join(characters))
        push!(latex, replace(join(inputs, ", "), "\\\\", "\\"))
        push!(desc, join(unicode_names, " + "))
    end
    return code, unicode, latex, desc
end

open("$filename.txt","w") do f
    print(f, """
        $filename.txt  LaTeX-to-Unicode reference table

        ===================================================================
        LATEX-TO-UNICODE REFERENCE TABLE    *L2U-ref* *julia-vim-L2U-reference*

          Note: This file is autogenerated from the script '$(basename(Base.source_path()))'
          The symbols are based on the documentation of Julia version $VERSION
          See |julia-vim| for the LaTeX-to-Unicode manual.

        """)

    col_headers = ["Code point(s)", "Character(s)", "Tab completion sequence(s)", "Unicode name(s)"]

    code, unicode, latex, desc =
        table_entries(
            tab_completions(
                Base.REPLCompletions.latex_symbols,
                Base.REPLCompletions.emoji_symbols
                ),
            unicode_data()
            )

    cw = max(length(col_headers[1]), maximum(map(length, code)))
    uw = max(length(col_headers[2]), maximum(map(length, unicode)))
    lw = max(length(col_headers[3]), maximum(map(length, latex)))
    dw = max(length(col_headers[4]), maximum(map(length, desc)))

    print_padded(c, u, l, d) = println(f, rpad(c, cw), " ", rpad(u, uw), " ", rpad(l, lw), " ", d)

    print_padded(col_headers[1:3]..., col_headers[4] * "~")
    print_padded("-"^cw, "-"^uw, "-"^lw, "-"^dw)

    for (c, u, l, d) in zip(code, unicode, latex, desc)
        print_padded(c, u, l, d)
    end
    print_padded("-"^cw, "-"^uw, "-"^lw, "-"^dw)

    println(f)
    println(f, "  vim", # separated to avoid an error from vim which otherwise tries to parse this line
               ":tw=$(cw+uw+lw+dw+3):et:ft=help:norl:")
end

max_stars_repo_name: myokoym/dotfiles

max_issues_repo_name: myokoym/dotfiles

max_forks_repo_name: myokoym/dotfiles


text:
# Univariate Options
function optimize(f::F,
     lower::T,
     upper::T;
     method = Brent(),
     rel_tol::Real = sqrt(eps(T)),
     abs_tol::Real = eps(T),
     iterations::Integer = 1_000,
     store_trace::Bool = false,
     show_trace::Bool = false,
     callback = nothing,
     show_every = 1,
     extended_trace::Bool = false) where {F<:Function, T <: AbstractFloat}
    show_every = show_every > 0 ? show_every : 1
    if extended_trace && callback == nothing
        show_trace = true
    end

    show_trace && print_header(method)

    optimize(f, lower, upper, method;
             rel_tol = T(rel_tol),
             abs_tol = T(abs_tol),
             iterations = iterations,
             store_trace = store_trace,
             show_trace = show_trace,
             show_every = show_every,
             callback = callback,
             extended_trace = extended_trace)
end

function optimize(f::F,
    lower::Union{Integer, Real},
    upper::Union{Integer, Real};
    kwargs...) where F<:Function

    T = promote_type(typeof(lower/1), typeof(upper/1))
    optimize(f,
             T(lower),
             T(upper);
             kwargs...)
end

function optimize(f::F,
    lower::Union{Integer, Real},
    upper::Union{Integer, Real},
    method::Union{Brent, GoldenSection};
    kwargs...) where F<:Function

    T = promote_type(typeof(lower/1), typeof(upper/1))
    optimize(f,
             T(lower),
             T(upper),
             method;
             kwargs...)
end


text:
export rhotrace,
    rhotrace!,
    steptrace,
    steptrace!,
    radiustrace,
    radiustrace!,
    objtrace,
    objtrace!,
    gradtrace,
    gradtrace!

function log_display(v, name)
    if (m = minimum(v)) > 0.0
        return v
    elseif m ≈ 0.0
        @warn "$name: shifted by $(eps(Float64)) to display on Log Scale"
        return v .+ eps(Float64)
    end
    @warn "$name: shifted by $m to display on Log Scale"
    return v .+ m
end

"""
    objtrace

Provide any number of terminal simulations as arguments and `objtrace` displays
a plot of each arguments' objective function values at each succussful iterate.

To add a simulation to the current `objtrace` plot use `objtrace!`.
"""
function objtrace end

"""
    gradtrace

Provide any number of terminal simulations as arguments and `gradtrace` displays
a plot of each arguments' normed gradient values at each succussful iterate.

To add a simulation to the current `gradtrace` plot use `gradtrace!`.
"""
function gradtrace end

"""
    radiustrace

Provide any number of terminal simulations as arguments and `radiustrace` displays
a plot of each arguments' trust region subproblem radius values at each succussful iteration.

To add a simulation to the current `radiustrace` plot use `radiustrace!`.
"""
function radiustrace end

"""
    steptrace

Provide any number of terminal simulations as arguments and `steptrace` displays
a plot of each arguments' step measure at each succussful iteration.

To add a simulation to the current `steptrace` plot use `steptrace!`.
"""
function steptrace end

"""
    rhotrace

Provide any number of terminal simulations as arguments and `rhotrace` displays
a plot of each arguments' ρ ratio at each successful iteration.

To add a simulation to the current `rhotrace` plot use `rhotrace!`.
"""
function rhotrace end

@userplot ObjTrace

@recipe function f(result::ObjTrace)
    num_args = length(result.args)

    m = model(result.args[1])

    y_lab = (isa(formula(m), Missing) || num_args > 1 ? "\$ f \$" : formula(m))

    xguide := "Iterates (gHS evaluations)"

    yguide := y_lab

    yguidefontsize := 7

    xguidefontsize := 7

    title --> "Objective Trace"

    titlefontsize --> 11

    yaxis --> :log

    legend --> :bottomright

    @series begin
        x = []

        labels = Matrix{String}(undef, 1, num_args)

        # Current idea is to accept trace, backend, simulation objects (may need to include model, driver in simulation.jl)
        for i ∈ 1:num_args

            s = result.args[i]

            qn_update = string(QN_update(driver(s)))

            s_update = "6.1." * string(S_update(driver(s)))[end]

            two_w = samples(driver(s))

            dim = dimension(model(s))

            id = name(model(s))

            v = log_display(f_vals(s), id)

            labels[
                1,
                i,
            ] = "\$ f \\quad \\mathrm{$id}: ~ n=$dim \\quad w=$(Int(two_w/2.0)) \\quad (\\mathrm{$s_update}) \\quad \\mathrm{$qn_update} \$"

            push!(x, v)
        end

        label --> labels

        x
    end
end

@userplot GradTrace

@recipe function f(result::GradTrace)
    num_args = length(result.args)

    m = model(result.args[1])

    xguide := "Iterates (gHS evaluations)"

    yguide := "\$ ||\\nabla f ||\$"

    yguidefontsize := 7

    xguidefontsize := 7

    title --> "Normed Gradient Trace"

    titlefontsize --> 11

    yaxis --> :log

    legend --> :bottomright
    @series begin
        x = []

        labels = Matrix{String}(undef, 1, num_args)

        for i ∈ 1:num_args

            s = result.args[i]

            qn_update = string(QN_update(driver(s)))

            s_update = "6.1." * string(S_update(driver(s)))[end]

            two_w = samples(driver(s))

            dim = dimension(model(s))

            id = name(model(s))

            v = log_display(∇f_norms(s), id)

            labels[
                1,
                i,
            ] = "\$ ||\\nabla f|| ~ \\mathrm{$id}: ~ n=$dim \\quad w=$(Int(two_w/2.0)) \\quad (\\mathrm{$s_update}) \\quad \\mathrm{$qn_update} \$"

            push!(x, v)
        end

        label --> labels

        x
    end
end

@userplot RadiusTrace

@recipe function f(result::RadiusTrace)

    xguide := "Successful Steps"

    yguide := "\$ \\Delta_k \$"

    yguidefontsize := 7

    xguidefontsize := 7

    title --> "Trust-Region Subproblem Radius Trace"

    titlefontsize --> 11

    yaxis --> :log

    legend --> :bottomright

    @series begin
        x = []

        labels = Matrix{String}(undef, 1, length(result.args))

        for i ∈ 1:length(result.args)

            s = result.args[i]

            qn_update = string(QN_update(driver(s)))

            s_update = "6.1." * string(S_update(driver(s)))[end]

            two_w = samples(driver(s))

            dim = dimension(model(s))

            id = name(model(s))

            v = log_display(Δ_vals(s), id)

            labels[
                1,
                i,
            ] = "\$ \\Delta_k ~ \\mathrm{$id}: ~ n=$dim \\quad w=$(Int(two_w/2.0)) \\quad (\\mathrm{$s_update}) \\quad \\mathrm{$qn_update} \$"

            push!(x, v)
        end

        label --> labels

        x
    end
end

@userplot StepTrace

@recipe function f(result::StepTrace)

    xguide := "Successful Steps"

    yguide := "\$ ||p_k|| \$"

    yguidefontsize := 7

    xguidefontsize := 7

    title --> "Step Distance Trace"

    titlefontsize --> 11

    yaxis --> :log

    legend --> :bottomright

    @series begin
        x = []

        labels = Matrix{String}(undef, 1, length(result.args))

        for i ∈ 1:length(result.args)

            s = result.args[i]

            qn_update = string(QN_update(driver(s)))

            s_update = "6.1." * string(S_update(driver(s)))[end]

            two_w = samples(driver(s))

            dim = dimension(model(s))

            id = name(model(s))

            v = log_display(p_norms(s), id)

            labels[
                1,
                i,
            ] = "\$ ||p_k|| ~ \\mathrm{$id}: ~ n=$dim \\quad w=$(Int(two_w/2.0)) \\quad (\\mathrm{$s_update}) \\quad \\mathrm{$qn_update} \$"

            push!(x, v)
        end

        label --> labels

        x
    end

end

@userplot RhoTrace

@recipe function f(result::RhoTrace)

    xguide := "Successful Steps"

    yguide := "\$ \\rho \$"

    yguidefontsize := 7

    xguidefontsize := 7

    title --> "Actual Reduction to Model Reduction Ratio Trace"

    titlefontsize --> 11

    yaxis --> :log

    legend --> :bottomright

    @series begin
        x = []

        labels = Matrix{String}(undef, 1, length(result.args))

        for i ∈ 1:length(result.args)

            s = result.args[i]

            qn_update = string(QN_update(driver(s)))

            s_update = "6.1." * string(S_update(driver(s)))[end]

            two_w = samples(driver(s))

            dim = dimension(model(s))

            id = name(model(s))

            v = log_display(ρ_vals(s), id)

            labels[
                1,
                i,
            ] = "\$ \\rho ~ \\mathrm{$id}: ~ n=$dim \\quad w=$(Int(two_w/2.0)) \\quad (\\mathrm{$s_update}) \\quad \\mathrm{$qn_update} \$"

            push!(x, v)
        end

        label --> labels

        x
    end
end


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

using Markdown
using InteractiveUtils

# ╔═╡ 02bfa078-d62b-11ea-15df-d701431829b9
begin
	using Pkg
	Pkg.activate("..")
	Pkg.instantiate()
	using Unitful
	using Measurements, UnitfulAstro, Parameters, Markdown
	using PhysicalConstants.CODATA2018: G, k_B, m_u, σ
	const amu, k = m_u, k_B
end;

# ╔═╡ c9ac27ee-dac0-11ea-2a8c-2d144b034a82
md"""
# Exoplanet Calculator 🪐
"""

# ╔═╡ b2286b26-dac2-11ea-1ce0-c7da562aa641
md"""
This calculator uses host star and orbital parameters from the literature to self-consistently calculate derived values used for estimating the observed signal (ΔD) produced by its planet's atmosphere.

For proper rendering of unicode symbols, we recommend using either Chrome of Firefox to view this webpage. The source code for this page is available [here](https://github.com/icweaver/exocalc).
"""

# ╔═╡ 19b35ef4-dac3-11ea-2d25-97e5482ff6a0
md"### Literature values"

# ╔═╡ 07db65d6-dd99-11ea-103b-33d0317af127
md"Reference values from [NASA Exoplanet Archive](https://exoplanetarchive.ipac.caltech.edu/cgi-bin/DisplayOverview/nph-DisplayOverview?objname=HAT-P-23%20b)"

# ╔═╡ 7cff6dfc-dd9f-11ea-1fdf-7b6aaa9435b4
md"### Results"

# ╔═╡ a35f2936-e41d-11ea-3bdc-3b347410a558
md"Scroll to see more:"

# ╔═╡ 38a61304-e0fe-11ea-14b2-17d9b9e13c7b
md"### Calculate parameters"

# ╔═╡ 0b6821a4-dac3-11ea-27d7-911521f0d3c0
md"### How parameters were calculated"

# ╔═╡ f8281da6-dd9f-11ea-1b6c-d32702215397
md"""
The calculator first checks if there are missing input parameters and then calls the appropriate function to calculate them for each study. The resulting derived parameters are then calculated from them.

Everything was done with a "star first" approach, meaning that all stellar parameters were determined first, and then the planet parameters were determined self-consistently from that. If conflicting parameters are given, the calculator will try to give priority to using direct observables to perform calculations and error otherwise.

For transparency, the inputs used for each calculation are shown in parenthesis next to each parameter.
"""

# ╔═╡ 49f75dea-dda0-11ea-1a85-bbdd4750b878
md"""
These are the functions used to calculate each parameter based on the combination of inputs given. No `if` statements or default "`None`" keyword arguments needed thanks to Julia's multiple dispatch! For ease of reference, these functions also return a string of the inputs used.
"""

# ╔═╡ 3f79c516-da77-11ea-1f6b-d3e7191a95d8
begin
	# Star radius
	get_Rₛ(Lₛ, Tₛ) = (Lₛ / (4.0π*σ*Tₛ^4))^(1//2), "Rₛ(Lₛ, Tₛ)"

	# Star-to-planet radius ratio
	get_RₚRₛ(Rₚ, Rₛ) = Rₚ / Rₛ, "RₚRₛ(Rₚ, Rₛ)"

	# Semi-major axis / Star density
	function get_aRₛ(ρₛ::Unitful.Density, P::Unitful.Time)
		return ((G * P^2 * ρₛ)/(3.0π))^(1//3), "aRₛ(ρₛ, P)"
	end
	get_aRₛ(a::Unitful.Length, Rₛ::Unitful.Length) = a / Rₛ, "aRₛ(a, Rₛ)"

	# Semi-major axis
    get_a(aRₛ, Rₛ) = aRₛ * Rₛ, "a(aRₛ, Rₛ)"

	# Impact parameter
	get_b(i, aRₛ) = aRₛ * cos(i), "b(i, aRₛ)"

	# Star density
	function get_ρₛ(P::Unitful.Time, aRₛ::Measurement)
        return (3.0π / (G * P^2)) * aRₛ^3, "ρₛ(P, aRₛ)"
	end
	function get_ρₛ(Mₛ::Unitful.Mass, Rₛ::Unitful.Length)
		return Mₛ / ((4.0/3.0)π * Rₛ^3), "ρₛ(Mₛ, Rₛ)"
	end

	# Star mass
	get_Mₛ(ρₛ, Rₛ) = ρₛ * (4.0/3.0) * π * Rₛ^3.0, "Mₛ(ρₛ, Rₛ)"

	# Star luminosity
	get_Lₛ(Tₛ, Rₛ) = 4.0π * Rₛ^2 * σ * Tₛ^4, "Lₛ(Tₛ, Rₛ)"

	# Star temperature
	get_Tₛ(Lₛ, Rₛ) = (L / (4.0π * Rₛ^2 * σ))^(1//4), "Tₛ(Lₛ, Rₛ)"

	# Planet mass
	get_Mₚ(K, i, P, Mₛ) =
		(K/sin(i)) * (P / (2.0π*G))^(1//3) * Mₛ^(2//3), "Mₚ(K, i, P, Mₛ)"

	# Planet radius
	get_Rₚ(RₚRₛ, Rₛ) =  RₚRₛ * Rₛ, "Rₚ(RₚRₛ, Rₛ)"

	# Planet density
	get_ρₚ(Mₚ, Rₚ) = Mₚ / ((4.0/3.0)π * Rₚ^3), "ρₚ(Mₚ, Rₚ)"

	# Star surface gravity
	get_gₛ(Mₛ, Rₛ) = G * Mₛ / Rₛ^2, "gₛ(Mₛ, Rₛ)"

	# Planet surface gravity
	get_gₚ(Mₚ, RₚRₛ, Rₛ) = G * Mₚ / (RₚRₛ^2 * Rₛ^2), "gₚ(Mₚ, RₚRₛ, Rₛ)"

	# Planet equilibrium temperature
	get_Tₚ(Tₛ, aRₛ, α) = Tₛ * (1.0 - α)^(1//4) * (0.5/aRₛ)^(1//2), "Tₚ(Tₛ, aRₛ, α)"

	# Planet scale height
	get_H(μ, Tₚ, gₚ) = k * Tₚ / (μ * gₚ), "H(μ, Tₚ, gₚ)"

	# Estimated signal from planet atmosphere
	get_ΔD(H, RₚRₛ, Rₛ) = 2.0 * H * RₚRₛ/Rₛ
end;

# ╔═╡ 8e5811ae-dd9e-11ea-127e-b9812511492b
md"### Input parameters from study"

# ╔═╡ db28dbd2-db12-11ea-28e4-2b6cf30bd102
@with_kw_noshow struct Study @deftype Union{Nothing, Quantity, Measurement}
	# Reference name (e.g. Ciceri et al. 2015)
	name::String = "Custom"

	# Star params
	Tₛ = nothing
	ρₛ = nothing
	Mₛ = nothing
	Rₛ = nothing
	gₛ = nothing
	Lₛ = nothing

	# Orbital params
	RₚRₛ = nothing
	aRₛ = nothing
	a = nothing
	b = nothing
	P = nothing
	K = nothing
	i = nothing

	# Planet params
	μ = nothing
	α = nothing
	Tₚ = nothing
	ρₚ = nothing
	Mₚ = nothing
	Rₚ = nothing
	gₚ = nothing
	N_scales::Float64 = 5.0 # Number of scale heights (default 5)
end;

# ╔═╡ 17302d74-d63b-11ea-3de3-49f0df0554ca
# Input params from studies to explore
studies = [
	Study(
		name = "WASP-43/b: Weaver et al. (2020)",
		μ	 = 2.0*amu,
		α	 = 0.0 ± 0.0,
		i	 = (1.433 ± 0.1)u"rad",
		P	 = (0.813473978 ± 3.5e-8)u"d",
		Tₛ	 = (4520 ± 120)u"K",
		Rₛ	 = (0.667 ± 0.010)u"Rsun",
		aRₛ  = 4.872 ± 0.14,
		Mₛ	 = (0.717 ± 0.025)u"Msun",
		Tₚ	 = (1440.0 ± 40.0)u"K",
		Mₚ	 = (2.052 ± 0.053)u"Mjup",
		K    = (551.7 ± 4.7)u"m/s",
		Rₚ	 =	(1.036 ± 0.012)u"Rjup",
	),
	Study(
		name = "HAT-P-23/b: Ciceri et al. (2015)",
		μ	 = 2.0*amu,
		α	 = 0.0 ± 0.0,
		K	 = (368.5 ± 17.6)u"m/s",
		i	 = (85.74 ± 0.95)u"°",
		P	 = (1.21288287 ± 0.00000017)u"d",
		RₚRₛ = 0.11616 ± 0.00081,
		Tₛ	 = (5885.0 ± 72.0)u"K",
		Rₛ	 = (1.089 ± 0.028)u"Rsun",
		aRₛ  = (4.5459 ± 0.0919),
	),
	Study(
		name = "HAT-P-23/b: Sada & Ramón-Fox (2016)",
		μ	 = 2.0*amu,
		α	 = 0.0 ± 0.0,
		K	 = (346.0 ± 21.0)u"m/s", # latest RV data, from B17
		i	 = (85.1 ± 1.5)u"°",  # latest RV data, from B17
		P	 = (1.212880 ± 0.000002)u"d", # latest transit data: (S&R16)
		RₚRₛ = 0.1113 ± 0.0010, # latest transit data: (S&R16)
		Tₛ	 = (5905.0 ± 80.0)u"K",
		Rₛ	 = (0.960±0.200)u"Rsun",
		aRₛ  = 4.26 ± 0.14,
	),
	Study(
		name = "HAT-P-23/b: Stassun et al. (2017, GAIA DR1)",
		μ	 = 2.0*amu,
		α	 = 0.0 ± 0.0,
		K	 = (368.5 ± 17.6)u"m/s",
		i	 = (85.1 ± 1.5)u"°",
		P	 = (1.212880 ± 0.000002)u"d",
		RₚRₛ = 0.1113 ± 0.0010,
		Tₛ	 = (5905.0 ± 80.0)u"K",
		ρₛ	 = (0.92 ± 0.18)u"g/cm^3",
		Rₛ	 = (0.960 ± 0.200)u"Rsun",
	),
	Study(
		name = "HAT-P-23/b: GAIA DR2",
		μ	 = 2.0*amu,
		α	 = 0.0 ± 0.0,
		K	 = (346.0 ± 21.0)u"m/s", # latest RV data, from B17
		i	 = (85.1 ± 1.5)u"°",  # latest RV data, from B17
		P	 = (1.2128867 ± 0.0000002)u"d", # latest transit data: (S&R16)
		RₚRₛ = 0.1113 ± 0.0010, # latest transit data: (S&R16)
		Tₛ	 = (5734.0 ± 100.0)u"K",
		Rₛ	 = (1.1858169 ± 0.0424133)u"Rsun",
		Lₛ	 = (10.0^(0.13656067 ± 0.00864667))u"Lsun",
		#Mₚ	 = (1.34 ± 0.59)u"Mjup", # DR1 mass, gives inconsistent ΔD=550ppm result
		aRₛ  = (4.26 ± 0.14), # latest transit data: (S&R16)
	),
	Study(
		name = "HAT-P-23/b: TICv8",
		μ	 = 2.0*amu,
		α	 = 0.0 ± 0.0,
		K	 = (346.0 ± 21.0)u"m/s", # latest RV data, from B17
		i	 = (85.1 ± 1.5)u"°",  # latest RV data, from B17
		P	 = (1.2128867 ± 0.0000002)u"d", # latest transit data: (S&R16)
		RₚRₛ = 0.1113 ± 0.0010, # latest transit data: (S&R16)
		Tₛ	 = (5918.230 ± 136.811)u"K",
		Rₛ	 = (1.1517600 ± 0.0596583)u"Rsun",
		ρₛ   = (0.99471000 ± 0.23240140)u"g/cm^3",
		Mₛ   = (1.078000 ± 0.136618)u"Msun",
		Lₛ   = (10.0^(0.1661873 ± 0.0191600))u"Lsun",
		gₛ   = (10.0^(4.3479600 ± 0.0819789))u"cm/s^2",
	),
	Study(
		name = "HAT-P-23/b: Set 1 (Gaia DR1)",

		# Weaver et al. (2020)
		i	 = (83.6 ± 0.3)u"°",
		P	 = (1.21289 ± 3.73975e-8)u"d",
		RₚRₛ = 0.11390 ± 0.0010,
		#ρₛ = (1.016 ± 0.0286)u"g/cm^3",

		# Literature/Fixed parameter inputs
		μ	 = 2.0*amu,
		α	 = 0.0 ± 0.0,
		K	 = (346.0 ± 21.0)u"m/s", # latest RV data, from B17
		Tₛ	 = (5905.0 ± 80.0)u"K", # GAIA DR1
		Rₛ	 = (0.960 ± 0.200)u"Rsun", # GAIA DR1
		ρₛ	 = (0.92 ± 0.18)u"g/cm^3", # GAIA DR1
	),
	Study(
		name = "HAT-P-23/b: Set 2 (TICv8)",

		# Weaver et al. (2020)
		i	 = (83.6 ± 0.3)u"°",
		P	 = (1.21289 ± 3.73975e-8)u"d",
		RₚRₛ = 0.11390 ± 0.0010,
		#ρₛ   = (1.016 ± 0.0286)u"g/cm^3",

		# Literature/Fixed parameter inputs
		μ	 = 2.0*amu,
		α	 = 0.0 ± 0.0,
		K	 = (346.0 ± 21.0)u"m/s", # latest RV data, from B17
		Tₛ	 = (5918.230 ± 136.811)u"K", # TICv8 (DR2)
		Rₛ	 = (1.1517600 ± 0.0596583)u"Rsun", # TICv8 (DR2)
		ρₛ   = (0.99471000 ± 0.23240140)u"g/cm^3", # TICv8 (DR2)

	),
	Study(
		name = "WASP-50b: Chakrabarty & Sengupta (2019)",
		i	 = (84.88 ± 0.27)u"°",
		P	 = (1.955100 ± 0.000005)u"d",
		RₚRₛ = 0.1390 ± 0.0006,
		μ	 = 2.0*amu,
		α	 = 0.0 ± 0.0,
		K	 = (256.6 ± 4.4)u"m/s",
		Tₛ	 = (5400 ± 100)u"K",
		Rₛ	 = (0.843 ± 0.031)u"Rsun",
		aRₛ  = 7.51 ± 0.10,
	),
	Study(
		name = "WASP-50b: C&S 2019 + TICv8",
		i	 = (84.88 ± 0.27)u"°",
		P	 = (1.955100 ± 0.000005)u"d",
		RₚRₛ = 0.1390 ± 0.0006,
		μ	 = 2.0*amu,
		α	 = 0.0 ± 0.0,
		K	 = (256.6 ± 4.4)u"m/s",
		Tₛ	 = (5518.0000 ± 128.3340)u"K", # TICv8
		Rₛ	 = (0.8734180 ± 0.0453186)u"Rsun", # TICv8
		ρₛ   = (2.0503245 ± 0.4676309)u"g/cm^3", # TICv8 (not self consistent)
		#Mₛ   = (0.96900000 ± 0.14486300)u"Msun", # TICv8
	),
	Study(
		name = "WASP-50b: TESS multi-sector + TICv8",
		i	 = (84.5157 ± 0.207)u"°",
		P	 = (1.9550931258 ± 3.963e-7)u"d",
		RₚRₛ = 0.1370858217 ± 0.0013727689,
		μ	 = 2.0*amu,
		α	 = 0.0 ± 0.0,
		K	 = (256.6 ± 4.4)u"m/s", # C&S 2019
		Tₛ	 = (5518.0000 ± 128.3340)u"K", # TICv8
		Rₛ	 = (0.8734180 ± 0.0453186)u"Rsun", # TICv8
		ρₛ   = (2.0503245 ± 0.4676309)u"g/cm^3", # TICv8
		#aRₛ  = 7.3624634217 ± 0.10479750132
	),
];

# ╔═╡ a8df7ad0-dd9e-11ea-2a6a-f16683371016
md"### Derived parameters"

# ╔═╡ bd752a9e-dd80-11ea-141c-779c5135d4d8
@with_kw_noshow struct Derived @deftype Tuple{Union{Nothing, Quantity, Measurement}, String}
	# Reference name (e.g. Ciceri et al. 2015)
	name::String = "Custom"

	# Star Params
	ρₛ
	gₛ
	Mₛ
	Rₛ
	Tₛ
	Lₛ

	#Orbital params
	RₚRₛ::Tuple{Measurement, String}
	P
	aRₛ::Tuple{Measurement, String}
	a
	b::Tuple{Measurement, String}
	K
	i

	# Planet params
	μ
	α::Tuple{Measurement, String}
	gₚ
	Mₚ
	Rₚ
	ρₚ
	Tₚ
	H

	# Signal
	N_scales::Float64
	ΔD::Measurement
end;

# ╔═╡ c01eb856-e0f9-11ea-01d5-07593189ce46
function calculate_params(st::Study)
	# Rₛ, Tₛ, Lₛ
	if !isnothing(st.Rₛ)
		Rₛ, inputs_Rₛ = st.Rₛ, "Rₛ(Rₛ)"
		if any((!isnothing).([st.Tₛ, st.Lₛ]))
			if all((!isnothing).([st.Tₛ, st.Lₛ]))
				Lₛ, inputs_Lₛ = st.Lₛ, "Lₛ(Lₛ)"
				Tₛ, inputs_Tₛ = st.Tₛ, "Tₛ(Tₛ)"
			elseif isnothing(st.Tₛ)
				Lₛ, inputs_Lₛ = st.Lₛ, "Lₛ(Lₛ)"
				Tₛ, inputs_Tₛ = get_Tₛ(Lₛ, Rₛ)
			else
				Tₛ, inputs_Tₛ = st.Tₛ, "Tₛ(Tₛ)"
				Lₛ, inputs_Lₛ = get_Lₛ(Tₛ, Rₛ)
			end
		else
			error("Rₛ was given. Lₛ or Tₛ must also be given.")
		end
	else
		if all((!isnothing).([st.Lₛ, st.Tₛ]))
			Lₛ, inputs_Lₛ = st.Lₛ, "Lₛ(Lₛ)"
			Tₛ, inputs_Tₛ = st.Tₛ, "Tₛ(Tₛ)"
			Rₛ, inputs_Rₛ = get_Rₛ(Lₛ, Tₛ)
		else
			error("Rₛ was not given. Lₛ and Tₛ must be given then.")
		end
	end

	# RₚRₛ and Rₚ
	if !isnothing(st.Rₚ)
		Rₚ, inputs_Rₚ = st.Rₚ, "Rₚ(Rₚ)"
		RₚRₛ, inputs_RₚRₛ = get_RₚRₛ(Rₚ, Rₛ)
	elseif !isnothing(st.RₚRₛ)
		RₚRₛ, inputs_RₚRₛ = st.RₚRₛ, "RₚRₛ(RₚRₛ)"
		Rₚ, inputs_Rₚ = get_Rₚ(RₚRₛ, Rₛ)
	else
		error("Please specify either RₚRₛ or Rₚ.")
	end

	# P
	if !isnothing(st.P)
		P, inputs_P = st.P, "P(P)"
	else
		error("Please specify a period.")
	end

	# ρₛ, aRₛ, a
	if all((!isnothing).([st.aRₛ, st.ρₛ]))
		error("Conflicting inputs. Only aRₛ or ρₛ can be given.")
	end
	if all((!isnothing).([st.aRₛ, st.b]))
		error("Conflicting inputs. Only aRₛ or b can be given.")
	end
	if !isnothing(st.ρₛ)
		ρₛ, inputs_ρₛ = st.ρₛ, "ρₛ(ρₛ)"
		aRₛ, inputs_aRₛ = get_aRₛ(ρₛ, P)
		a, inputs_a = get_a(aRₛ, Rₛ)
	elseif !isnothing(st.aRₛ)
		aRₛ, inputs_aRₛ = st.aRₛ, "aRₛ(aRₛ)"
		a, inputs_a = get_a(aRₛ, Rₛ)
		ρₛ, inputs_ρₛ = get_ρₛ(P, aRₛ)
	elseif !isnothing(st.a)
		a, inputs_a = st.a, "a(a)"
		aRₛ, inputs_aRₛ = get_aRₛ(a, Rₛ)
		ρₛ, inputs_ρₛ = get_ρₛ(P, aRₛ)
	# else
	# 	error("ρₛ or (aRₛ or a, or Mₛ) must be given for $(st.name)")
	end

	# Mₛ
	if !isnothing(st.Mₛ)
		Mₛ, inputs_Mₛ = st.Mₛ, "Mₛ(Mₛ)"
		ρₛ, inputs_ρₛ = get_ρₛ(Mₛ, Rₛ)
		aRₛ, inputs_aRₛ = get_aRₛ(ρₛ, P)
		a, inputs_a = get_a(aRₛ, Rₛ)
	else
		Mₛ, inputs_Mₛ = get_Mₛ(ρₛ, Rₛ)
	end

	# Calculate remaining params if not given/calculated
	i, inputs_i   =
	!isnothing(st.i) ? (st.i, "i(i)") : error("Must provide inclination (i).")
	K, inputs_K   =
	!isnothing(st.K) ? (st.K, "K(K)") : error("Must provide RV semi-amplitude (K).")
	α, inputs_α   =
	!isnothing(st.α) ? (st.α, "α(α)") : error("Must provide albedo (α).")
	b,  inputs_b  = !isnothing(st.b)  ? (st.b, "b(b)")   : get_b(i, aRₛ)
	Mₚ, inputs_Mₚ = !isnothing(st.Mₚ) ? (st.Mₚ, "Mₚ(Mₚ)") : get_Mₚ(K, i, P, Mₛ)
	Tₚ, inputs_Tₚ = !isnothing(st.Tₚ) ? (st.Tₚ, "Tₚ(Tₚ)") : get_Tₚ(Tₛ, aRₛ, α)
    gₛ, inputs_gₛ = !isnothing(st.gₛ) ? (st.gₛ, "gₛ(gₛ)")   : get_gₛ(Mₛ, Rₛ)
	gₚ, inputs_gₚ = !isnothing(st.gₚ) ? (st.gₚ, "gₚ(gₚ)") : get_gₚ(Mₚ, RₚRₛ, Rₛ)
	ρₚ, inputs_ρₚ = !isnothing(st.ρₚ) ? (st.ρₚ, "ρₚ(ρₚ)") : get_ρₚ(Mₚ, Rₚ)
	μ, inputs_μ =
	!isnothing(st.μ) ? (st.μ, "μ(μ)") :
	error("Must provide mean molecula weight (μ).")

	# Calculate signal
	H, inputs_H  = get_H(μ, Tₚ, gₚ)
	ΔD = get_ΔD(H, RₚRₛ, Rₛ)

	# Store results
	return Derived(
		# Study
		name = st.name,

		# Star Params
		ρₛ = (ρₛ, inputs_ρₛ),
		gₛ = (gₛ, inputs_gₛ),
		Mₛ = (Mₛ, inputs_Mₛ),
		Rₛ = (Rₛ, inputs_Rₛ),
		Tₛ = (Tₛ, inputs_Tₛ),
		Lₛ = (Lₛ, inputs_Lₛ),

		#Orbital params
		RₚRₛ = (RₚRₛ, inputs_RₚRₛ),
		P = (P, inputs_P),
		aRₛ = (aRₛ, inputs_aRₛ),
		a = (a, inputs_a),
		b = (b, inputs_b),
		K = (K, inputs_K),
		i = (i, inputs_i),

		# Planet params
		μ = (μ, inputs_μ),
		α = (α, inputs_α),
		gₚ = (gₚ, inputs_gₚ),
		Mₚ = (Mₚ, inputs_Mₚ),
		Rₚ = (Rₚ, inputs_Rₚ),
		ρₚ = (ρₚ, inputs_ρₚ),
		Tₚ = (Tₚ, inputs_Tₚ),
		H = (H, inputs_H),

		# Signal
		N_scales = st.N_scales,
		ΔD = ΔD,
	)
end;

# ╔═╡ a5a123d6-c0e7-4542-a231-6c50970d1470
const results = Derived[calculate_params(st) for st in studies]

# ╔═╡ c5c5ea28-dd9e-11ea-1f89-5b1371831177
md"### Display final results"

# ╔═╡ 33fc58d0-dbd9-11ea-3c45-83f4b5a2a818
function display_summary(d::Derived)
	md"""
	###### **$(d.name):**
	**Star Params** \
	$(d.Rₛ[2]) = $(d.Rₛ[1] |> u"Rsun") \
	$(d.Mₛ[2]) = $(d.Mₛ[1] |> u"Msun") \
	$(d.Tₛ[2]) = $(d.Tₛ[1] |> u"K") \
	$(d.Lₛ[2]) = $(d.Lₛ[1] |> u"Lsun") \
	$(d.ρₛ[2]) = $(d.ρₛ[1] |> u"g/cm^3") \
	$(d.gₛ[2])(cm/s²) = $(log10(ustrip(d.gₛ[1] |> u"cm/s^2")))

	**Orbital params** \
	$(d.K[2]) = $(d.K[1] |> u"m/s") \
	$(d.i[2]) = $(d.i[1] |> u"°") \
	$(d.RₚRₛ[2]) = $(d.RₚRₛ[1] |> NoUnits) \
	$(d.aRₛ[2]) = $(d.aRₛ[1] |> NoUnits) \
	$(d.P[2]) = $(d.P[1] |> u"d") \
	$(d.b[2]) = $(d.b[1] |> NoUnits)

	**Planet params** \
	$(d.μ[2]) = $(d.μ[1] |> u"u") \
	$(d.α[2]) = $(d.α[1]) \
	$(d.Rₚ[2]) = $(d.Rₚ[1] |> u"Rjup") \
	$(d.Mₚ[2]) = $(d.Mₚ[1] |> u"Mjup") \
	$(d.ρₚ[2]) = $(d.ρₚ[1] |> u"g/cm^3") \
	$(d.Tₚ[2]) = $(d.Tₚ[1] |> u"K") \
	$(d.gₚ[2]) = $(d.gₚ[1] |> u"m/s^2") \
	log10 g (cm/s^2) = $(d.gₚ[1] |> u"cm/s^2" |> ustrip |> log10) \
	$(d.H[2]) = $(d.H[1] |> u"km")

	**Signal at $(d.N_scales) scale heights** \
	ΔD = $(d.N_scales * uconvert(NoUnits, d.ΔD) * 1e6) ppm
	"""
end;

# ╔═╡ 6836afe2-6846-4582-8afd-c3d3a1e832ba
display_summary.(results)

# ╔═╡ 7db94ad6-dda1-11ea-2f33-1da144f1b7ad
md"Libraries for using things like physical constants and units."

# ╔═╡ Cell order:
# ╟─c9ac27ee-dac0-11ea-2a8c-2d144b034a82
# ╟─b2286b26-dac2-11ea-1ce0-c7da562aa641
# ╟─19b35ef4-dac3-11ea-2d25-97e5482ff6a0
# ╟─07db65d6-dd99-11ea-103b-33d0317af127
# ╠═17302d74-d63b-11ea-3de3-49f0df0554ca
# ╟─7cff6dfc-dd9f-11ea-1fdf-7b6aaa9435b4
# ╟─a35f2936-e41d-11ea-3bdc-3b347410a558
# ╠═6836afe2-6846-4582-8afd-c3d3a1e832ba
# ╟─38a61304-e0fe-11ea-14b2-17d9b9e13c7b
# ╠═a5a123d6-c0e7-4542-a231-6c50970d1470
# ╟─0b6821a4-dac3-11ea-27d7-911521f0d3c0
# ╟─f8281da6-dd9f-11ea-1b6c-d32702215397
# ╠═c01eb856-e0f9-11ea-01d5-07593189ce46
# ╟─49f75dea-dda0-11ea-1a85-bbdd4750b878
# ╠═3f79c516-da77-11ea-1f6b-d3e7191a95d8
# ╟─8e5811ae-dd9e-11ea-127e-b9812511492b
# ╠═db28dbd2-db12-11ea-28e4-2b6cf30bd102
# ╟─a8df7ad0-dd9e-11ea-2a6a-f16683371016
# ╠═bd752a9e-dd80-11ea-141c-779c5135d4d8
# ╟─c5c5ea28-dd9e-11ea-1f89-5b1371831177
# ╠═33fc58d0-dbd9-11ea-3c45-83f4b5a2a818
# ╟─7db94ad6-dda1-11ea-2f33-1da144f1b7ad
# ╠═02bfa078-d62b-11ea-15df-d701431829b9

max_stars_repo_name: icweaver/exocalc

max_issues_repo_name: icweaver/exocalc

max_forks_repo_name: icweaver/exocalc


text:
using Meganet
using Base.Test

for TYPE=[Float64,Float32]
    @testset "Tikhonov $TYPE" begin
    S = getTikhonovReg(TYPE)

    x = randn(TYPE,10)
    Sc,dS,d2S = regularizer(S,x)

    @test norm(x-dS)/norm(x) < 1e3*eps(TYPE)
    end
end


text:
using ViewerGL; GL = ViewerGL
using LinearAlgebraicRepresentation
Lar = LinearAlgebraicRepresentation

#	1D complex in 2D generated by Lar and PLaSM operators
circle = Lar.circle()()
GL.VIEW([ GL.GLFrame2, GL.GLGrid( circle...,GL.COLORS[1],1 ) ]);

intervals = Lar.qn(1)([1.5,1.5])
GL.VIEW([ GL.GLFrame2, GL.GLGrid( intervals...,GL.COLORS[1],1 ) ]);

#	2D complex in 3D generated by Lar and PLaSM operators
cylinders = Lar.larModelProduct([ circle, intervals ])
GL.VIEW([ GL.GLFrame2, GL.GLGrid( cylinders...,GL.COLORS[1],1 ) ]);

function cylinder(r,h,n,k)
	circle = Lar.circle(r)([n])
	interval = Lar.qn(k)([h])
	V,FV = Lar.larModelProduct([ circle, interval ])
	push!(FV, [k for k=1:2*n if k%2==0])
	push!(FV, [k for k=1:2*n if k%2==1])

	circle1 = circle[1],[[k] for k=1:size(circle[1],2)]
	interval0 = interval[1], [[k] for k=1:size(interval[1],2)]
	_,EV0 = Lar.larModelProduct([ circle, interval0 ])
	_,EV1 = Lar.larModelProduct([ circle1, interval ])
	EV = append!(EV0,EV1)
	return V,FV,EV
end

rod = Lar.apply(Lar.t(0,0,-1.5), cylinder(1,3,12,1))
triple = Lar.Struct([ rod, Lar.r(pi/2,0,0), rod, Lar.r(0,pi/2,0), rod ])
V,(_,EV,FV,_) = Lar.apply( Lar.s(2.5,2.5,2.5)*Lar.t(-.50,-.5,-.5), Lar.cuboidGrid([1,1,1],true) )
cube = V,FV,EV
W, EW, FW = catmullclark(V, EV, FV, 2)
sphere = Lar.apply( Lar.s( 1.4, 1.4, 1.4), (W,FW,EW) )
object = Lar.Struct([ triple, cube, sphere ])
V,FV,EV = Lar.struct2lar(object)

GL.VIEW([ GL.GLFrame, GL.GLLines(V,EV) ]);

cop_EV = Lar.coboundary_0(EV::Lar.Cells);
cop_EW = convert(Lar.ChainOp, cop_EV);
cop_FE = Lar.coboundary_1(V, FV::Lar.Cells, EV::Lar.Cells);
W = convert(Lar.Points, V');

V, copEV, copFE, copCF = Lar.Arrangement.spatial_arrangement( W, cop_EW, cop_FE)


text:
using GridFunctions
using StaticArrays
using Test

@testset "GridArray D=$D" for D in 1:3
    dom = Domain{Float64,D}(fill(0, SVector{D,Float64}), fill(1, SVector{D,Float64}))
    bb = BBox{D}(fill(0, SVector{D,Int}), fill(10, SVector{D,Int}))
    cc = fill(false, SVector{D,Bool})

    ga = GridArray{Float64,Nothing,D}(dom, bb, cc, fill(nothing, Tuple(size(bb))))

    test_AbstractGridComponent(ga)
end