text:
#    Copyright (c) 2021 J.A. Duffek
#    Copyright (c) 2000 D.M. Spink
#
#    This program is free software; you can redistribute it and/or modify
#    it under the terms of the GNU General Public License as published by
#    the Free Software Foundation; either version 2 of the License, or
#    (at your option) any later version.
#
#    This program is distributed in the hope that it will be useful,
#    but WITHOUT ANY WARRANTY; without even the implied warranty of
#    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
#    GNU General Public License for more details.
#
#    You should have received a copy of the GNU General Public License
#    along with this program; if not, see <http://www.gnu.org/licenses/>.

"""
    vecnorm_toolbox(vec::Union{Matrix{F}, Array{F,3}}
                   )::Union{Matrix{F}, Array{F,3}} where {F<:AbstractFloat}

or

    vecnorm_toolbox(vec::Vector{F})::Vector{F} where {F<:AbstractFloat}

Normalises the array of vectors, returning the unit vectors.

# Arguments:
- `vec`: An array of column vectors represented by a matrix of size `(dim,nv)`,
  where `dim` is the dimension of the vector and `nv` the number of vectors.

# Output:
- `nvec`: Normalised vectors, matrix the same size as `vec`.

# Examples:
```julia
julia> nvec = vecnorm_toolbox(vec)
```

Normalise the two vectors `[0.0;2.0;1.3]` and `[1.5;3.4;2.3]`.
```julia
julia> nvec = vecnorm_toolbox([0.0 1.5; 2.0 3.4; 1.3 2.3])
```
"""
function vecnorm_toolbox(vec::Union{Matrix{F}, Array{F,3}}
                        )::Union{Matrix{F},Array{F,3}} where {F<:AbstractFloat}
return vec ./ sqrt.(sum(vec.^2,dims=1));
end # vecnorm_toolbox

# TODO: look if this is even necessary for this toolbox (using LinearAlgebra)
# overloaded function for vector input to get a float instead of a vector
function vecnorm_toolbox(vec::Vector{F})::Vector{F} where {F<:AbstractFloat}
return vec ./ sqrt(vec' * vec);
end # vecnorm_toolbox


text:
module Combinatorics

using ...Prerequisites.TypeTraits: corder,indtosub

export AbstractCombinatorics
export Combinations,DulCombinations,Permutations,DulPermutations

"""
    AbstractCombinatorics{M,C}

Abstract combinatoric algorithms.
"""
abstract type AbstractCombinatorics{M,C} end
Base.eltype(::Type{<:AbstractCombinatorics{M,C}}) where {M,C}=NTuple{M,C|>eltype}
@generated gettuple(contents,inds,::Val{M}) where M=Expr(:tuple,[:(contents[inds[$i]]) for i=1:M]...)

"""
    Combinations{M}(contents::C) where {M,C}

Combinations of M elements from contents. Duplicates are not allowed.
"""
struct Combinations{M,C} <: AbstractCombinatorics{M,C}
    contents::C
    N::Int
    Combinations{M}(contents::C) where {M,C}=new{M,C}(contents,contents|>length)
end
Base.length(c::Combinations{M}) where M=binomial(c.N,M)
Base.iterate(c::Combinations{M}) where M=M>c.N ? nothing : M==0 ? ((),[c.N+2]) : (gettuple(c.contents,1:M,Val(M)),nextmstate!(collect(1:M),c.N,M))
Base.iterate(c::Combinations{M},state) where M=state[1]>c.N-M+1 ? nothing : (gettuple(c.contents,state,Val(M)),nextmstate!(state,c.N,M))
function nextmstate!(state::Vector{Int},N::Int,M::Int)
    for i=M:-1:1
        state[i]+=1
        state[i]>(N-(M-i)) && continue
        for j=i+1:M
            state[j]=state[j-1]+1
        end
        break
    end
    state
end

"""
    DulCombinations{M}(contents::C) where {M,C}

Combinations of M elements from contents. Duplicates are allowed.
"""
struct DulCombinations{M,C} <: AbstractCombinatorics{M,C}
    contents::C
    N::Int
    DulCombinations{M}(contents::C) where {M,C}=new{M,C}(contents,contents|>length)
end
Base.length(c::DulCombinations{M}) where M=binomial(c.N+M-1,c.N-1)
Base.iterate(c::DulCombinations{M}) where M=M==0 ? ((),[c.N+1]) : (gettuple(c.contents,ntuple(i->1,M),Val(M)),nextdmstate!(collect(ntuple(i->1,M)),c.N,M))
Base.iterate(c::DulCombinations{M},state) where M=state[1]>c.N ? nothing : (gettuple(c.contents,state,Val(M)),nextdmstate!(state,c.N,M))
function nextdmstate!(state::Vector{Int},N::Int,M::Int)
    for i=M:-1:1
        state[i]+=1
        state[i]>N && continue
        for j=i+1:M
            state[j]=state[j-1]
        end
        break
    end
    state
end

"""
    Permutations{M}(contents::C) where {M,C}

Permutations of M elements from contents. Duplicates are allowed.
"""
struct Permutations{M,C} <: AbstractCombinatorics{M,C}
    contents::C
    N::Int
    Permutations{M}(contents::C) where {M,C}=new{M,C}(contents,contents|>length)
end
Base.length(p::Permutations{M}) where M=(0<=M<=p.N) ? prod((p.N-M+1):p.N) : 0
Base.iterate(p::Permutations{M}) where M=(p.N==0 && M>0 || 0<p.N<M) ? nothing : (state=collect(1:p.N);(gettuple(p.contents,state,Val(M)),nextpstate!(state,p.N,M)))
Base.iterate(p::Permutations{M},state) where M=(p.N==0 && M>0 || 0<p.N<max(state[1],M)) ? nothing : (gettuple(p.contents,state,Val(M)),nextpstate!(state,p.N,M))
function nextpstate!(state::Vector{Int},N::Int,M::Int)
    M<=0 && return [N+1]
    M<N && (j=M+1; while j<=N && state[M]>=state[j]; j+=1; end)
    if M<N && j<=N
        state[M],state[j]=state[j],state[M]
    else
        M<N && reverse!(state,M+1)
        i=M-1
        while i>=1 && state[i]>=state[i+1];i-=1; end
        if i>0
            j=N
            while j>i && state[i]>=state[j];j-=1; end
            state[i],state[j]=state[j],state[i]
            reverse!(state,i+1)
        else
            state[1]=N+1
        end
    end
    return state
end

"""
    DulPermutations{M}(contents::C) where {M,C}

Permutations of M elements from contents. Duplicates are not allowed.
"""
struct DulPermutations{M,C} <: AbstractCombinatorics{M,C}
    contents::C
    N::Int
    DulPermutations{M}(contents::C) where {M,C}=new{M,C}(contents,contents|>length)
end
Base.length(p::DulPermutations{M}) where M=p.N^M
Base.iterate(p::DulPermutations{M},state=1) where M=state>length(p) ? nothing : (indtosub(ntuple(i->p.N,M),state,corder),state+1)

end # module

ext: jl


text:
##

using ACE
using Printf, Test, LinearAlgebra, StaticArrays
using ACE: evaluate, evaluate_d, Rn1pBasis, Ylm1pBasis,
      PositionState, Product1pBasis, State, ACEConfig,
      SymmetricBasis
using Random: shuffle
using ACEbase.Testing: dirfdtest, fdtest, print_tf, test_fio, println_slim
using ACE.OrthPolys: transformed_jacobi

##

@info("======= Testing XScal1pBasis ======= ")

maxdeg = 5
trans = ACE.Transforms.IdTransform()
P = transformed_jacobi(2*maxdeg, trans, 1.0, 0.0; pin = 0, pcut = 0)

Bsel = ACE.SimpleSparseBasis(3, maxdeg)

##

B1p = ACE.xscal1pbasis(P, (k = 1:maxdeg, m = 0:maxdeg), :u; label = "Xkm")
ACE.init1pspec!(B1p, Bsel)
ACE.fill_rand_coeffs!(B1p, randn)

##

@info("check symbols")
println_slim(@test sort(ACE.symbols(B1p)) == sort([:k, :m]) )

@info("check `get_index`")
for ntest = 1:30
   local idx, b
   idx = rand(1:length(B1p))
   b = B1p.spec[idx]
   print_tf(@test ACE.get_index(B1p, b) == idx)
end
println()

@info("check fio")
println_slim(@test all( test_fio(B1p; warntype=false) ))

##

@info("evaluation test ")
X = State(u = rand())
B = evaluate(B1p, X)
dB = evaluate_d(B1p, X)
println_slim(@test ACE.evaluate_ed(B1p, X) == (B, dB) )

##

@info("Finite-difference tests at a few random points")
for ntest = 1:10
   local c, F, dF
   c = rand(length(B1p))
   F = u -> ACE.contract(c, evaluate(B1p, State(u = u)))
   dF = u -> ACE.contract(c, getproperty.(evaluate_d(B1p, State(u = u)), :u))
   print_tf(@test all( fdtest(F, dF, rand() - 0.5; verbose=false) ))
end
println()

##

@info("check diag coeffs")
ACE.fill_diag_coeffs!(B1p, :k)
u = rand()
J = evaluate(B1p.P, u)
B = evaluate(B1p, State(u=u))
println_slim(@test all(B[ib] == J[b.k] for (ib, b) in enumerate(B1p.spec)))

##

@info("incorporate into product basis")

Xka_u = ACE.xscal1pbasis(P, (k = 1:maxdeg, a = 0:maxdeg),:u; label = "Xka_u")
ACE.init1pspec!(Xka_u, Bsel)
Pa_v = ACE.Scal1pBasis(:v, nothing, :a, P, "Pa_v")
B1p = Xka_u * Pa_v
ACE.init1pspec!(B1p, Bsel)
println_slim(@test B1p["Xka_u"] === Xka_u)
println_slim(@test B1p["Pa_v"] === Pa_v)

rand_uv_state() = State(u = rand(), v = rand())
X = rand_uv_state()

evaluate(B1p, X)

@show length(B1p)
@show length(Xka_u)
@show length(Pa_v)
@info("sparsify")
ACE.sparsify!(B1p, ACE.get_spec(B1p))
@show length(B1p)
@show length(Xka_u)
@show length(Pa_v)

##

@info("XScal with norm(rr) as input")
B1p = ACE.xscal1pbasis(P, (k = 1:maxdeg, m = 0:maxdeg), ACE.GetNorm{:rr}();
                       label = "Xkm")
ACE.init1pspec!(B1p, Bsel)
ACE.fill_rand_coeffs!(B1p, randn)

@info("check get_val, get_val_d")
X = State( rr = ACE.rand_sphere() * (0.5+rand()/2) )
println_slim(@test ACE.getval(X, B1p) ≈ norm(X.rr))
println_slim(@test ACE.getval_d(X, B1p).rr ≈ X.rr/norm(X.rr))

@info("check get_val, get_val_d")
B = evaluate(B1p, X)
dB = evaluate_d(B1p, X)
println_slim(@test ACE.evaluate_ed(B1p, X) == (B, dB))

##

@info("Finite-difference tests at a few random points")
for ntest = 1:10
   local c, F, dF
   c = randn(length(B1p))
   rr0 = Vector(ACE.rand_sphere() * (0.5 + rand()/0.5))
   F = rr -> ACE.contract(c, evaluate(B1p, State(rr = SVector{3}(rr))))
   dF = rr -> ACE.contract(c, getproperty.(evaluate_d(B1p, State(rr = SVector{3}(rr))), :rr)) |> Vector
   dF(rr0)
   print_tf(@test all( fdtest(F, dF, rr0; verbose=false) ))
end
println()

##

@info("Compatibility Test")

maxdeg = 15
Bsel = ACE.SimpleSparseBasis(3, maxdeg)
trans = ACE.Transforms.IdTransform()
P = transformed_jacobi(maxdeg, trans, 1.0, 0.0; pin = 0, pcut = 0)

@info(" ... Scal1pBasis compatibility")
Bu = ACE.Scal1pBasis(:u, nothing, :k, P)
Bu_x = ACE.xscal1pbasis(P, (k = 1:maxdeg,), ACE.GetVal{:u}(), label = "Bu_x")
ACE.init1pspec!(Bu_x, Bsel)
Bu_x.coeffs[:,:] += I

for ntest = 1:30
   X = ACE.State(u = rand())
   print_tf(@test ACE.evaluate(Bu, X) ≈ ACE.evaluate(Bu_x, X) )
end
println()

@info(" ... Rn1pBasis compatibility")
Br = ACE.Rn1pBasis(P; label = "Rn", varsym = :rr, nsym = :n)
Br_x = ACE.xscal1pbasis(P, (n = 1:maxdeg,), ACE.GetNorm{:rr}(); label = "Rn_x")
ACE.init1pspec!(Br_x, Bsel)
Br_x.coeffs[:,:] += I

for ntest = 1:30
   rr = ACE.rand_sphere() * rand()
   X = ACE.State(rr = rr)
   print_tf(@test ACE.evaluate(Br, X) ≈ ACE.evaluate(Br_x, X) )
end
println()


text:
import Compat: readstring

let
    # Ensure that open cannot find the file "foo"
    @test !isfile("foo")
    @test_throws SystemError open("foo")

    # Note that patches always have higher precedence. The patch `open(::Any)`
    # will be preferred over the original `open(::AbstractString)` for `open("foo")`
    patch = @patch open(name) = IOBuffer("bar")
    apply(patch) do
        @test readstring(@mock open("foo")) == "bar"

        # The `open(::Any)` patch could result in unintended (or intended) consequences
        @test readstring(@mock open(`echo helloworld`)) == "bar"
    end

    # Better to be specific with your patches
    patch = @patch open(name::AbstractString) = IOBuffer("bar")
    apply(patch) do
        @test readstring(@mock open("foo")) == "bar"

        # The more specific `open(::AbstractString)` patches only a single method
        io, pobj = (@mock open(`echo helloworld`))
        @test readstring(io) == "helloworld\n"
    end
end


text:
using MICS
using DataFrames
using CSV
using Base.Test

β = 1.6773985789

states = Vector{State}()
push!( states, State( readtable("log_1.dat",separator=' '), @prop(β*:E1), @prop(β*:E2) ) )
push!( states, State( readtable("log_2.dat",separator=' '), @prop(β*:E2) ) )
push!( states, State( readtable("log_3.dat",separator=' '), @prop(β*:E3) ) )
push!( states, State( readtable("log_4.dat",separator=' '), @prop(β*:E4) ) )

@time mixture = Mixture( states; title = "Test", verbose=true )

#@time f, δf = freeEnergies( mixture )


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

using Markdown
using InteractiveUtils

# This Pluto notebook uses @bind for interactivity. When running this notebook outside of Pluto, the following 'mock version' of @bind gives bound variables a default value (instead of an error).
macro bind(def, element)
    quote
        local el = $(esc(element))
        global $(esc(def)) = Core.applicable(Base.get, el) ? Base.get(el) : missing
        el
    end
end

# ╔═╡ bc2fdc0c-04b2-11eb-26dd-ab2d62794853
# Include package.
begin
	using Pkg
	Pkg.activate(joinpath(Pkg.envdir(), "dev-env"))
	using FractalTools
	using Plots
	plotly()
	default(:label, "")
	default(:lw, 2)
end

# ╔═╡ c5132b30-04b2-11eb-173c-898862565f53
# Define function
begin
	f(x) = 1 / 2 - (x - 1 / 2)^2 + 1
	g(x) = x^2
end

# ╔═╡ 924bbe12-04b3-11eb-3ab4-f1d84faeb86f
begin
	xi, dx, xf = 0, 0.25,  1
	x = collect(xi : dx : xf)
	y = f.(x)
end

# ╔═╡ ee35dc64-04b3-11eb-22c3-5fc92810c3f1

@bind d html"<input type=range min=0 max=0.2 step=0.01>"

# ╔═╡ 30386406-04b4-11eb-2e1b-2d0eee87f5ed
@bind h  html"<input type=range min=0 max=0.2 step=0.01>"

# ╔═╡ 321a63d2-04b4-11eb-283d-9fa41a2534ef
@bind l  html"<input type=range min=0 max=0.2 step=0.01>"

# ╔═╡ 47922bdc-04b4-11eb-239b-411078a8b27b
@bind m  html"<input type=range min=0 max=0.2 step=0.01>"

# ╔═╡ 5f9872c4-04b4-11eb-162f-9f998adb0c9f
@bind niter html"<input type=range min=1 max=30 step=1>"

# ╔═╡ 508bc310-04b4-11eb-3d4c-e78390b8b8a1
# Construct interpolants
begin
	z = g.(x)
	dv =  d * ones(length(x) - 1)
	hv =  h * ones(length(x) - 1)
	lv =  l * ones(length(x) - 1)
	mv =  m * ones(length(x) - 1)
	interp = interpolate(x, y, d = dv, niter=niter)
	hiddeninterp = hiddenfif(x, y, z, d = dv, h = hv, l = lv, m = mv, niter=niter)
end

# ╔═╡ 34596672-04b5-11eb-1f41-f781db5f981f
d, h, l, m, niter

# ╔═╡ 99eb66ee-04b3-11eb-34e3-5f2bcf5f383f
begin
	xt = collect(xi : 0.1dx : xf)
	plt = plot()
	scatter!(x, y, marker=(:circle, 3), label="data")
	plot!(xt, f.(xt), label="f")
	plot!(xt, interp.(xt), label="finterp")
	plot!(xt, getindex.(hiddeninterp.(xt), 1), label="finterphidden")
end

# ╔═╡ Cell order:
# ╠═bc2fdc0c-04b2-11eb-26dd-ab2d62794853
# ╠═c5132b30-04b2-11eb-173c-898862565f53
# ╠═924bbe12-04b3-11eb-3ab4-f1d84faeb86f
# ╠═508bc310-04b4-11eb-3d4c-e78390b8b8a1
# ╠═ee35dc64-04b3-11eb-22c3-5fc92810c3f1
# ╠═30386406-04b4-11eb-2e1b-2d0eee87f5ed
# ╠═321a63d2-04b4-11eb-283d-9fa41a2534ef
# ╠═47922bdc-04b4-11eb-239b-411078a8b27b
# ╠═5f9872c4-04b4-11eb-162f-9f998adb0c9f
# ╠═34596672-04b5-11eb-1f41-f781db5f981f
# ╠═99eb66ee-04b3-11eb-34e3-5f2bcf5f383f


text:
###############################################################################
#
# Update canopy matrices
#
###############################################################################
"""
    canopy_matrices!(
                leaves::Array{LeafBios{FT},1},
                can_opt::CanopyOpticals{FT}
    ) where {FT<:AbstractFloat}

Compute scattering coefficient matrices for direct and diffuse light given
    geometry dependent overall extinction coefficients and pigment dependent
    leaf reflectance and transmission (computed via fluspect). This function
    has to be called before [`short_wave!`](@ref) can be used.
- `leaves` Array of [`LeafBios`](@ref) type struct
- `can_opt` [`CanopyOpticals`](@ref) type struct
"""
function canopy_matrices!(
            leaves::Array{LeafBios{FT},1},
            can_opt::CanopyOpticals{FT}
) where {FT<:AbstractFloat}
    # 1. unpack values
    @unpack ddb, ddf, dob, dof, sdb, sdf, sob, sof = can_opt;

    # 2. Calculation of reflectance
    nLayer = size(can_opt.sigb)[2];
    @inbounds for i=1:nLayer
        if length(leaves)>1
            τ_SW = leaves[i].τ_SW;
            ρ_SW = leaves[i].ρ_SW;
        else
            τ_SW = leaves[1].τ_SW;
            ρ_SW = leaves[1].ρ_SW;
        end

        #CF: Right now, canopy geometry is the same everywhere,
        #    can be easily extended to layers as well.
        @inbounds for j=1:size(can_opt.sigb, 1)
            # diffuse backscatter coefficient for diffuse incidence
            can_opt.sigb[j,i] = ddb * ρ_SW[j] + ddf * τ_SW[j];
            # diffuse forwardscatter coefficient for diffuse incidence
            can_opt.sigf[j,i] = ddf * ρ_SW[j] + ddb * τ_SW[j];
            # diffuse backscatter coefficient for specular incidence
            can_opt.sb[j,i]   = sdb * ρ_SW[j] + sdf * τ_SW[j];
            # diffuse forwardscatter coefficient for specular incidence
            can_opt.sf[j,i]   = sdf * ρ_SW[j] + sdb * τ_SW[j];
            # directional backscatter  coefficient for diffuse incidence
            can_opt.vb[j,i]   = dob * ρ_SW[j] + dof * τ_SW[j];
            # directional forwardscatter coefficient for diffuse incidence
            can_opt.vf[j,i]   = dof * ρ_SW[j] + dob * τ_SW[j];
            # bidirectional scattering coefficent (directional-directional)
            can_opt.w[j,i]    = sob * ρ_SW[j] + sof * τ_SW[j];
        end
    end

    # 3. attenuation
    can_opt.a .= 1 .- can_opt.sigf;

    return nothing
end

max_stars_repo_name: bbuman/Land

max_issues_repo_name: bbuman/Land

max_forks_repo_name: bbuman/Land


text:
import Base.+, Base.-, Base.*, Base./, Base.convert, Base.zero, Base.one

export SampledSpectrum, SingleLine, NoLight, nolight, make_CIE_table

using Pkg
using DelimitedFiles

struct SpectrumStyle <: Broadcast.AbstractArrayStyle{1} end
struct LineStyle <: Broadcast.BroadcastStyle end

# ------------------------------------------------
# Load CIE spectra

const CIE_PATH = joinpath(dirname(pathof(GRAYTR)), "..", "data", "lin2012xyz2e_1_7sf.csv")

const raw_CIE_data = readdlm(CIE_PATH, ',')
const N_CIE = size(raw_CIE_data)[1]

const CIE_DATA = SMatrix{4, N_CIE, Float64}(raw_CIE_data')
const CIE_YINT = sum(CIE_DATA[3,:])

struct CIE_Table
    data::Array{Float64,2}
    yint::Float64
end

CIE_Table(data) = CIE_Table(data, sum(data[:,2]))

# ------------------------------------------------
# Utility functions used when integrating spectral data

function find_bin(data, x)
    N_DATA = size(data)[1]
    for i = 1:N_DATA-1
        if data[i,1] <= x && data[i+1,1] > x
            return i
        end
    end
    nothing
end

function get_midpoint(DATA, x, n)
    a = DATA[n,1]
    b = DATA[n+1,1]
    t = (x - a) / (b - a)
    y = lerp(t, DATA[n,2], DATA[n+1,2])
    return t,y
end

function integrate(X, Y)
    deltas = X[2:end] - X[1:end-1]
    means = (Y[1:end-1] + Y[2:end]) / 2
    return sum(deltas .* means)
end

"""
    compute_bin(DATA; low, high ; averaged = false)

Compute the bin value by integrating over a (hopefully) higher-resolution spectrum.

- `DATA` is a `(N, 2)` array with wavelengths (in nm) in the first column and values in the second.
- `low` and `high` are the wavelength limits of the bin.
- `averaged` optionally returns an average over the spectrum instead of an integral
  (i.e. it divides the result by the bin length).

"""
function compute_bin(DATA, low, high ; averaged = false)
    nlow = find_bin(DATA, low)
    nhigh = find_bin(DATA, high)

    if nlow == nothing || nhigh == nothing
        return 0.0
    end

    if nlow == nhigh
        # There is only one data value in the range.
        # Return it if average, or return it times the interval if non-averaged
        low_t, low_y = get_midpoint(DATA, low, nlow)
        high_t, high_y = get_midpoint(DATA, high, nhigh)
        return 0.5 * (low_y + high_y) * (high - low)
    end

    low_t, low_y = get_midpoint(DATA, low, nlow)
    high_t, high_y = get_midpoint(DATA, high, nhigh)

    low_part =  0.5 * (low_y + DATA[nlow+1,2]) * (DATA[nlow+1,1] - low)
    high_part = 0.5 * (high_y + DATA[nlow,2]) * (high - DATA[nhigh,1])

    mid_part = integrate(DATA[nlow+1 : nhigh, 1], DATA[nlow+1 : nhigh, 2])

    S = low_part + mid_part + high_part

    return averaged ? S / (high - low) : S
end

# ------------------------------------------------
# The NoLight dummy spectrum

struct NoLight <: Spectrum end

const nolight = NoLight()

@inline isblack(N::NoLight) = true
@inline to_XYZ(S::NoLight) = (0.0, 0.0, 0.0)
to_XYZ(S::NoLight, C::CIE_Table) = (0.0, 0.0, 0.0)

Base.broadcastable(N::NoLight) = N

@inline interpolate(N::NoLight, x) = nolight
@inline integrate(N::NoLight) = 0.0

@inline Base.broadcasted(::typeof(+), s::NoLight, t::NoLight) = NoLight()
@inline Base.broadcasted(::typeof(*), n::NoLight, t::NoLight) = NoLight()

@inline Base.broadcasted(::typeof(*), n::Number, t::NoLight) = NoLight()
@inline Base.broadcasted(::typeof(*), t::NoLight, n::Number) = NoLight()
@inline Base.broadcasted(::typeof(+), n::Number, t::NoLight) = n
@inline Base.broadcasted(::typeof(+), t::NoLight, n::Number) = n

@inline Base.broadcasted(::typeof(*), t::NoLight, s::Spectrum) = NoLight()
@inline Base.broadcasted(::typeof(*), s::Spectrum, t::NoLight) = NoLight()
@inline Base.broadcasted(::typeof(+), t::NoLight, s::Spectrum) = s
@inline Base.broadcasted(::typeof(+), s::Spectrum, t::NoLight) = s

@inline Base.broadcasted(::typeof(*), t::NoLight, s::Broadcast.Broadcasted{SpectrumStyle,N}) where N = NoLight()
@inline Base.broadcasted(::typeof(+), t::NoLight, s::Broadcast.Broadcasted{SpectrumStyle,N}) where N = s

# ------------------------------------------------
# The SampledSpectrum type

struct SampledSpectrum <: Spectrum
    low::Int64
    high::Int64
    values::Array{Float64,1}
    delta::Float64
end

"""
    SampledSpectrum(low::Int64, high::Int64, values::Vector{Float64})

Create a new `SampledSpectrum` object with a wavelength range from `low` to `high`
nanometres and binned spectrum values given by `values`.

"""
function SampledSpectrum(low::Int64, high::Int64, values::Vector{Float64})
    delta = (high-low) / length(values)
    SampledSpectrum(low, high, values, delta)
end

Base.maximum(S::SampledSpectrum) = Base.maximum(S.values)
Base.minimum(S::SampledSpectrum) = Base.minimum(S.values)

isblack(S::SampledSpectrum) = Base.maximum(S.values) ≈ 0.0

# Various functions
Base.length(S::SampledSpectrum) = Base.length(S.values)
Base.size(S::SampledSpectrum) = Base.size(S.values)
Base.getindex(S::SampledSpectrum, inds::Vararg{Int,N}) where N = S.values[inds...]
Base.getindex(S::SampledSpectrum, inds::CartesianIndex{1}) = S.values[inds]
Base.setindex!(S::SampledSpectrum, val, inds::Vararg{Int,N}) where N = S.values[inds...] = val
Base.setindex!(S::SampledSpectrum, val, inds::CartesianIndex{1}) = S.values[inds] = val

Base.collect(S::SampledSpectrum) = S.values

Base.broadcastable(S::SampledSpectrum) = S

# Broadcasting style for SampledSpectrum
Base.BroadcastStyle(::Type{<:SampledSpectrum}) = SpectrumStyle()
Base.BroadcastStyle(::SpectrumStyle, ::SpectrumStyle) = SpectrumStyle()
Base.BroadcastStyle(::SpectrumStyle, ::Broadcast.AbstractArrayStyle{0}) = SpectrumStyle()
Base.BroadcastStyle(::SpectrumStyle, ::Broadcast.DefaultArrayStyle{0}) = SpectrumStyle()

function Base.similar(bc::Broadcast.Broadcasted{SpectrumStyle}, ::Type{ElType}) where ElType
    flat = Broadcast.flatten(bc)
    first = find_first(flat)
    l = first.low
    h = first.high
    N = length(first.values)

    for x in flat.args
        if isa(x, SampledSpectrum)
            if (x.low != l) || (x.high != h) || (length(x.values) != N)
                error("Spectral shape mismatch!")
            end
        end
    end

    return SampledSpectrum(first.low, first.high, similar(Array{ElType}, axes(bc)))
end

find_first(bc::Base.Broadcast.Broadcasted) = find_first(bc.args)
find_first(args::Tuple) = find_first(find_first(args[1]), Base.tail(args))
find_first(args::Tuple{}) = nothing
find_first(x) = x
find_first(a::SampledSpectrum, rest) = a
find_first(::Any, rest) = find_first(rest)

@inline isSS(S::Any) = false
@inline isSS(S::SampledSpectrum) = true

integrate(S::SampledSpectrum) = sum(S.values)

function interpolate(S::SampledSpectrum, a::Real)
    if a <= S.low || a >= S.high
        return 0.0
    end

    i = Int64(fld(a - S.low, S.delta))
    return S.values[i+1]
#    t = (a - i * S.delta) / S.delta
#    return @inbounds lerp(t, S.values[i+1], S.values[i+2])
end

function make_CIE_table(S::SampledSpectrum)
    N = length(S)
    data = zeros(N,3)
    delta = (S.high - S.low) / N
    CIE_lambda = raw_CIE_data[:,1]
    for n = 1:N
        # (low_x, high_x) are the wavelength limits of the bin; compute average over that
        low_x = S.low + (n-1) * delta
        high_x = S.low + n * delta
        data[n,1] = compute_bin(hcat(CIE_lambda, raw_CIE_data[:,2]), low_x, high_x ; averaged = true)
        data[n,2] = compute_bin(hcat(CIE_lambda, raw_CIE_data[:,3]), low_x, high_x ; averaged = true)
        data[n,3] = compute_bin(hcat(CIE_lambda, raw_CIE_data[:,4]), low_x, high_x ; averaged = true)
    end
    return CIE_Table(data)
end

function to_XYZ(S::SampledSpectrum, C::CIE_Table)
    x =y = z = 0.0
    @inbounds for i = 1:length(S)
        x += S.values[i] * C.data[i,1]
        y += S.values[i] * C.data[i,2]
        z += S.values[i] * C.data[i,3]
    end
    return x/C.yint, y/C.yint, z/C.yint
end

function to_XYZ(S::SampledSpectrum)
    x, y, z = 0.0, 0.0, 0.0
    @inbounds for i = 1:N_CIE
        lam = CIE_DATA[1,i]
        s = interpolate(S, lam)
        x += s * CIE_DATA[2,i]
        y += s * CIE_DATA[3,i]
        z += s * CIE_DATA[4,i]
    end
    return x/CIE_YINT, y/CIE_YINT, z/CIE_YINT
end

# ------------------------------------------------
# Single line spectra

"""
    SingleLine(wavelength, value)

Create a `SingleLine` spectrum representing a delta-distribution of a value at a certain
wavelength.

"""
struct SingleLine <: Spectrum
    wavelength::Float64
    value::Float64
end

isblack(S::SingleLine) = S.value ≈ 0.0

# Addition of SingleLine objects
function +(S1::SingleLine, S2::SingleLine)
    if S1.wavelength != S2.wavelength
        error("Spectrum wavelengths don't match:\n$S1 \n $S2")
    end
    return SingleLine(S1.wavelength, S1.value + S2.value)
end
+(S::SingleLine, c::Real) = SingleLine(S.wavelength, S.value + c)
+(c::Real, S::SingleLine) = S+c
-(c::Real, S::SingleLine) = SingleLine(S.wavelength, S.value - c)

# Multiplication of SingleLine objects
function *(S1::SingleLine, S2::SingleLine)
    if S1.wavelength != S2.wavelength
        SingleLine(S1.wavelength, 0.0)
    end
    return SingleLine(S1.wavelength, S1.value * S2.value)
end

Base.broadcastable(S::SingleLine) = S

# Multiplication broadcasting

@inline function Base.broadcasted(::typeof(*), S1::SingleLine, S2::SingleLine)
    if S1.wavelength != S2.wavelength
        return SingleLine(S1.wavelength, 0.0)
    end
    return SingleLine(S1.wavelength, S1.value * S2.value)
end

@inline function Base.broadcasted(::typeof(*), S1::SingleLine, S2::SampledSpectrum)
    return SingleLine(S1.wavelength, S1.value * interpolate(S2, S1.wavelength))
end
@inline function Base.broadcasted(::typeof(*), S2::SampledSpectrum, S1::SingleLine)
    return SingleLine(S1.wavelength, S1.value * interpolate(S2, S1.wavelength))
end
@inline function Base.broadcasted(::typeof(*), S1::SingleLine, n::Number)
    return SingleLine(S1.wavelength, S1.value * n)
end
@inline function Base.broadcasted(::typeof(*), n::Number, S1::SingleLine)
    return SingleLine(S1.wavelength, S1.value * n)
end

@inline function Base.broadcasted(::typeof(/), S1::SingleLine, n::Number)
    return SingleLine(S1.wavelength, S1.value / n)
end

@inline function Base.broadcasted(::typeof(+), S1::SingleLine, S2::SingleLine)
    if S1.wavelength != S2.wavelength
        return SingleLine(S1.wavelength, 0.0)
    end
    return SingleLine(S1.wavelength, S1.value + S2.value)
end

@inline function Base.broadcasted(::typeof(+), S1::SingleLine, S2::SampledSpectrum)
    return SingleLine(S1.wavelength, S1.value + interpolate(S2, S1.wavelength))
end
@inline function Base.broadcasted(::typeof(+), S2::SampledSpectrum, S1::SingleLine)
    return SingleLine(S1.wavelength, S1.value + interpolate(S2, S1.wavelength))
end
@inline function Base.broadcasted(::typeof(+), S1::SingleLine, n::Number)
    return SingleLine(S1.wavelength, S1.value + n)
end
@inline function Base.broadcasted(::typeof(+), n::Number, S1::SingleLine)
    return SingleLine(S1.wavelength, S1.value + n)
end

@inline Base.broadcasted(::typeof(*), S::SingleLine, B::Broadcast.Broadcasted{SpectrumStyle,N}) where N = S .* materialize(B)
@inline Base.broadcasted(::typeof(+), S::SingleLine, B::Broadcast.Broadcasted{SpectrumStyle,N}) where N = S .+ materialize(B)

function *(S1::SingleLine, S2::SampledSpectrum)
    return SingleLine(S1.wavelength, S1.value * interpolate(S2, S1.wavelength))
end
*(S1::SampledSpectrum, S2::SingleLine) = S2*S1

*(c::Real, S::SingleLine) = SingleLine(S.wavelength, c*S.value)
*(S::SingleLine, c::Real) = c*S

/(S::SingleLine, c::Real) = SingleLine(S.wavelength, S.value/c)

interpolate(S::SingleLine, wavelength::Real) = S.wavelength ≈ wavelength ? S.value : 0.0

integrate(S::SingleLine) = S.value

function to_XYZ(S::SingleLine)
    x,y,z = 0.0, 0.0, 0.0
    for i = 1:N_CIE-1
        if CIE_DATA[1,i] < S.wavelength && CIE_DATA[1,i+1] >= S.wavelength
            t = (CIE_DATA[1,i+1] - S.wavelength) / (CIE_DATA[1,i+1] - CIE_DATA[1,i])
            x += S.value * lerp(t, CIE_DATA[2,i], CIE_DATA[2,i+1])
            y += S.value * lerp(t, CIE_DATA[3,i], CIE_DATA[3,i+1])
            z += S.value * lerp(t, CIE_DATA[4,i], CIE_DATA[4,i+1])
            break
        end
    end
    return x/CIE_YINT, y/CIE_YINT, z/CIE_YINT
end

broadcast(f, S::SingleLine, y...) = SingleLine(S.wavelength, f(S.value, y...))

const Adobe_RGB = reshape([2.0413690, -0.5649464, -0.3446944,-0.9692660,  1.8760108,  0.0415560, 0.0134474, -0.1183897,  1.0154096], (3,3))'
const CIE_RGB = reshape([2.3706743, -0.9000405, -0.4706338,-0.5138850,  1.4253036,  0.0885814, 0.0052982, -0.0146949,  1.0093968], (3,3))'

# ------------------------------------------------
# Conversion between XYZ and RGB

XYZtoRGB(xyz::Vector) = CIE_RGB * xyz

function XYZtoRGB(x, y, z)
    return 2.3706743*x - 0.9000405*y - 0.4706338*z,
           -0.5138850*x + 1.4253036*y + 0.0885814*z,
            0.0052982*x - 0.0146949*y + 1.0093968*z
end

max_stars_repo_name: dronir/GRAYTR

max_issues_repo_name: dronir/GRAYTR

max_forks_repo_name: dronir/GRAYTR


text:
using OpenQuantumTools, OrdinaryDiffEq, Plots
H = DenseHamiltonian([(s)->1-s, (s)->s], -[σx, σz]/2)

# this plot recipe is for conveniently plotting the spectrum of the Hamiltonian
# the first 3 arguments are: the Hamiltonian, the grid `s` and the levels to keep
plot(H, 0:0.01:1, 2, linewidth=2)

H_h = DenseHamiltonian([(s)->1-s, (s)->s], -[σx, σz]/2, unit=:h)
H_ħ = DenseHamiltonian([(s)->1-s, (s)->s], -[σx, σz]/2, unit=:ħ)
println("Setting unit to :h")
@show evaluate(H_h, 0.5)
println("Setting unit to :ħ")
@show evaluate(H_ħ, 0.5);

H_h(0.5)

coupling = ConstantCouplings(["Z"])

η = 1e-4
fc = 4
T = 16
bath = Ohmic(η, fc, T)

info_freq(bath)

p1 = plot(bath, :γ, range(0,20,length=200), label="", size=(800, 400), linewidth=2)
p2 = plot(bath, :S, range(0,20,length=200), label="", size=(800, 400), linewidth=2)
plot(p1, p2, layout=(1,2), left_margin=3Plots.Measures.mm)

# Hamiltonian
H = DenseHamiltonian([(s)->1-s, (s)->s], -[σx, σz]/2, unit=:ħ)
# initial state
u0 = PauliVec[1][1]
# coupling
coupling = ConstantCouplings(["Z"], unit=:ħ)
# bath
bath = Ohmic(1e-4, 4, 16)
annealing = Annealing(H, u0; coupling=coupling, bath=bath)

tf = 10*sqrt(2)
@time sol = solve_schrodinger(annealing, tf, alg=Tsit5(), retol=1e-4)
# The following line of code is a convenient recipe to plot the instantaneous population during the evolution.
# It currently only supports Hamiltonians with an annealing parameter s = t/tf from 0 to 1.
# The third argument can be either a list or a number. When it is a list, it specifies the energy levels to plot (starting from 0); when it is a number, it specifies the total number of levels to plot.
plot(sol, H, [0], 0:0.01:tf, linewidth=2, xlabel = "t (ns)", ylabel="\$P_G(t)\$")

sol(0.5)

# You need to solve the unitary first before trying to solve Redfield equation
@time U = solve_unitary(annealing, tf, alg=Tsit5(), abstol=1e-8, retol=1e-8);
@time solve_von_neumann(annealing, tf, alg=Tsit5(), abstol=1e-8, retol=1e-8);

tf = 10*sqrt(2)
U = solve_unitary(annealing, tf, alg=Tsit5(), abstol=1e-8, retol=1e-8);
sol = solve_redfield(annealing, tf, U; alg=Tsit5(), abstol=1e-8, retol=1e-8)
plot(sol, H, [0], 0:0.01:tf, linewidth=2, xlabel="t (ns)", ylabel="\$P_G(t)\$")

tf = 10*sqrt(2)
@time sol = solve_ame(annealing, tf; alg=Tsit5(), ω_hint=range(-6, 6, length=100), reltol=1e-4)
plot(sol, H, [0], 0:0.01:tf, linewidth=2, xlabel="t (ns)", ylabel="\$P_G(t)\$")

tf = 5000
@time sol_ame = solve_ame(annealing, tf; alg=Tsit5(), ω_hint=range(-6, 6, length=100), reltol=1e-6)
plot(sol_ame, H, [0], 0:1:tf, linewidth=2, xlabel="t (ns)", ylabel="\$P_G(t)\$")

tf = 5000
# total number of trajectories
num_trajectories = 3000
# construct the `EnsembleProblem`
# `safetycopy` needs to be true because the current trajectories implementation is not thread-safe.
prob = build_ensembles(annealing, tf, :ame, ω_hint=range(-6, 6, length=200), safetycopy=true)
# to use multi-threads, you need to start the Julia kernel with multiple threads
# julia --threads 8
sol = solve(prob, Tsit5(), EnsembleThreads(), trajectories=num_trajectories, reltol=1e-6, saveat=range(0,tf,length=100))

t_axis = range(0,tf,length=100)
dataset = []
for t in t_axis
    w, v = eigen_decomp(H, t/tf)
    push!(dataset, [abs2(normalize(so(t))' * v[:, 1]) for so in sol])
end

# the following codes calculate the instantaneous ground state population and the corresponding error bars by averaging over all the trajectories:

pop_mean = []
pop_sem = []
for data in dataset
    p_mean = sum(data) / num_trajectories
    p_sem = sqrt(sum((x)->(x-p_mean)^2, data)) / num_trajectories
    push!(pop_mean, p_mean)
    push!(pop_sem, p_sem)
end

scatter(t_axis, pop_mean, marker=:d, yerror=2*pop_sem, label="Trajectory", markersize=6)
plot!(sol_ame, H, [0], t_axis, linewidth=2, label="Non-trajectory")
xlabel!("t (ns)")
ylabel!("\$P_G(s)\$")

using HOQSTTutorials
HOQSTTutorials.tutorial_footer(WEAVE_ARGS[:folder],WEAVE_ARGS[:file])


text:
import YAML
config = YAML.load(open("config.yaml"))

# Read the array from a text file

# 2D matrix of shape (npoints, nparams) spanning from [0,1]
nolh = readdlm("hypercube/array.dat")

# Map the ranges of parameters in config.yaml onto the Latin Hypercube

# Read emulator, and convert this to a 1D array of starting points, and
# a 1D array of widths

parnames = config["emulator"]["parameters"]

starting = Float64[config["emulator"]["ranges"][parname][1] for parname in parnames]
deltas = Float64[(config["emulator"]["ranges"][parname][2] - config["emulator"]["ranges"][parname][1]) for parname in parnames]

# Broadcast these to a proper 2D array
# Do this addition and broadcasting in column-major
shifts = (nolh' .* deltas)
parameters = shifts .+ starting

# Shift back to row-major
parameters = parameters'

# Now, write this parameter list to disk
writedlm("parameters.dat", parameters)

max_stars_repo_name: iancze/DiskEmulator

max_issues_repo_name: iancze/DiskEmulator

max_forks_repo_name: iancze/DiskEmulator


text:
abstract type LookupCanvas <: Canvas end

function lookup_encode end
function lookup_decode end

@inline pixel_width(c::LookupCanvas) = c.pixel_width
@inline pixel_height(c::LookupCanvas) = c.pixel_height
@inline origin_x(c::LookupCanvas) = c.origin_x
@inline origin_y(c::LookupCanvas) = c.origin_y
@inline width(c::LookupCanvas) = c.width
@inline height(c::LookupCanvas) = c.height
@inline grid(c::LookupCanvas) = c.grid
@inline colors(c::LookupCanvas) = c.colors
@inline nrows(c::LookupCanvas) = size(grid(c), 2)
@inline ncols(c::LookupCanvas) = size(grid(c), 1)

function CreateLookupCanvas(
        ::Type{T},
        char_width::Int,
        char_height::Int;
        origin_x::Number = 0.,
        origin_y::Number = 0.,
        width::Number = 1.,
        height::Number = 1.,
        min_char_height::Int = 5,
        min_char_width::Int = 2) where {T <: LookupCanvas}
    width  > 0 || throw(ArgumentError("width has to be positive"))
    height > 0 || throw(ArgumentError("height has to be positive"))
    char_width  = max(char_width, min_char_width)
    char_height = max(char_height, min_char_height)
    pixel_width  = char_width * x_pixel_per_char(T)
    pixel_height = char_height * y_pixel_per_char(T)
    grid   = fill(0x00, char_width, char_height)
    colors = fill(nothing, char_width, char_height)
    T(grid, colors, pixel_width, pixel_height,
      Float64(origin_x), Float64(origin_y),
      Float64(width), Float64(height))
end

function pixel!(c::T, pixel_x::Int, pixel_y::Int, color::UserColorType) where {T <: LookupCanvas}
    0 <= pixel_x <= pixel_width(c) || return c
    0 <= pixel_y <= pixel_height(c) || return c
    pixel_x = pixel_x < pixel_width(c) ? pixel_x : pixel_x - 1
    pixel_y = pixel_y < pixel_height(c) ? pixel_y : pixel_y - 1
    cw, ch = size(grid(c))
    tmp = pixel_x / pixel_width(c) * cw
    char_x = floor(Int, tmp) + 1
    char_x_off = (pixel_x % x_pixel_per_char(T)) + 1
    if char_x < round(Int, tmp, RoundNearestTiesUp) + 1 && char_x_off == 1
        char_x = char_x + 1
    end
    char_y = floor(Int, pixel_y / pixel_height(c) * ch) + 1
    char_y_off = (pixel_y % y_pixel_per_char(T)) + 1
    grid(c)[char_x, char_y] |= lookup_encode(c)[char_x_off, char_y_off]
    force = !(color isa Symbol)  # don't attempt to blend colors if they have been explicitly specified
    set_color!(c.colors, char_x, char_y, crayon_256_color(color); force=force)
    c
end

function printrow(io::IO, c::LookupCanvas, row::Int)
    0 < row <= nrows(c) || throw(ArgumentError("Argument row out of bounds: $row"))
    y = row
    for x in 1:ncols(c)
        print_color(colors(c)[x,y], io, lookup_decode(c)[grid(c)[x,y] + 1])
    end
end


text:
### This file contains convenience functions for constructing LoopSets.

function Base.copyto!(ls::LoopSet, q::Expr)
    q.head === :for || throw("Expression must be a for loop.")
    add_loop!(ls, q)
end

function add_ci_call!(q::Expr, f, args, syms, i)
    call = Expr(:call, f)
    for arg ∈ @view(args[2:end])
        if arg isa Core.SSAValue
            push!(call.args, syms[arg.id])
        else
            push!(call.args, arg)
        end
    end
    push!(q.args, Expr(:(=), syms[i], call))
end

function substitute_broadcast(q::Expr)
    ci = first(Meta.lower(LoopVectorization, q).args).code
    nargs = length(ci)-1
    ex = Expr(:block,)
    syms = [gensym() for _ ∈ 1:nargs]
    for n ∈ 1:nargs
        ciₙ = ci[n]
        ciₙargs = ciₙ.args
        f = first(ciₙargs)
        if ciₙ.head === :(=)
            push!(ex.args, Expr(:(=), f, syms[((ciₙargs[2])::Core.SSAValue).id]))
        elseif f === GlobalRef(Base, :materialize!)
            add_ci_call!(ex, lv(:vmaterialize!), ciₙargs, syms, n)
        elseif f === GlobalRef(Base, :materialize)
            add_ci_call!(ex, lv(:vmaterialize), ciₙargs, syms, n)
        else
            add_ci_call!(ex, f, ciₙargs, syms, n)
        end
    end
    ex
end

function LoopSet(q::Expr)
    q = SIMDPirates.contract_pass(q)
    ls = LoopSet()
    copyto!(ls, q)
    resize!(ls.loop_order, num_loops(ls))
    ls
end

"""
    @avx

Annotate a `for` loop, or a set of nested `for` loops whose bounds are constant across iterations, to optimize the computation. For example:

    function AmulBavx!(C, A, B)
        @avx for m ∈ 1:size(A,1), n ∈ 1:size(B,2)
            Cₘₙ = zero(eltype(C))
            for k ∈ 1:size(A,2)
                Cₘₙ += A[m,k] * B[k,n]
            end
            C[m,n] = Cₘₙ
        end
    end

The macro models the set of nested loops, and chooses a

It may also apply to broadcasts:

```jldoctest
julia> a = rand(100);

julia> b = @avx exp.(2 .* a);

julia> c = similar(b);

julia> @avx @. c = exp(2a);

julia> b ≈ c
true
```

"""
macro avx(q)
    q2 = if q.head === :for
        lower(LoopSet(q))
    else# assume broadcast
        substitute_broadcast(q)
    end
    esc(q2)
end


text:
module OneHotArrays

using Adapt
using ChainRulesCore
using GPUArrays
using LinearAlgebra
using MLUtils
using NNlib

export onehot, onehotbatch, onecold,
       OneHotArray, OneHotVector, OneHotMatrix, OneHotLike

include("array.jl")
include("onehot.jl")
include("linalg.jl")

end


text:
using BinaryBuilder

cuda_version = v"10.2.89"   # NOTE: could be less specific

sources_linux = [
    ArchiveSource("https://developer.nvidia.com/compute/machine-learning/cudnn/secure/7.6.5.32/Production/10.2_20191118/cudnn-10.2-linux-x64-v7.6.5.32.tgz",
                  "600267f2caaed2fd58eb214ba669d8ea35f396a7d19b94822e6b36f9f7088c20")
]
sources_windows = [
    ArchiveSource("https://developer.nvidia.com/compute/machine-learning/cudnn/secure/7.6.5.32/Production/10.2_20191118/cudnn-10.2-windows10-x64-v7.6.5.32.zip",
                  "fba812f60c61bc33b81db06cd55e8d769774d036186571d724295c71c9936064")
]

include("../common.jl")


text:
# Following the PowerModels convention exports everything except internal
# symbols, which are defined as those whose name starts with an underscore.
# If you don't want all of these symbols in your environment, then use
# `import` instead of `using`.

# Do not add items to this exclude list. Instead, rename them with an underscore.
const _EXCLUDE_SYMBOLS = [Symbol(@__MODULE__), :eval, :include]

for sym in names(@__MODULE__, all=true)
    sym_string = string(sym)
    if sym in _EXCLUDE_SYMBOLS || startswith(sym_string, "_") || startswith(sym_string, "@_")
        continue
    end
    if !(Base.isidentifier(sym) || (startswith(sym_string, "@") &&
         Base.isidentifier(sym_string[2:end])))
       continue
    end
    #println("$(sym)")
    @eval export $sym
end


text:
# ---
# title: 862. Shortest Subarray with Sum at Least K
# id: problem862
# author: Indigo
# date: 2021-02-17
# difficulty: Hard
# categories: Binary Search, Queue
# link: <https://leetcode.com/problems/shortest-subarray-with-sum-at-least-k/description/>
# hidden: true
# ---
#
# Return the **length** of the shortest, non-empty, contiguous subarray of `A`
# with sum at least `K`.
#
# If there is no non-empty subarray with sum at least `K`, return `-1`.
#
#
#
# **Example 1:**
#
#
#
#     Input: A = [1], K = 1
#     Output: 1
#
#
# **Example 2:**
#
#
#
#     Input: A = [1,2], K = 4
#     Output: -1
#
#
# **Example 3:**
#
#
#
#     Input: A = [2,-1,2], K = 3
#     Output: 3
#
#
#
#
# **Note:**
#
#   1. `1 <= A.length <= 50000`
#   2. `-10 ^ 5 <= A[i] <= 10 ^ 5`
#   3. `1 <= K <= 10 ^ 9`
#
#
## @lc code=start
using LeetCode

function shortest_subarray(A::Vector{Int}, K::Int)
    dq = Deque{Int}()
    prex = fill(0, length(A) + 1)
    cumsum!(@view(prex[2:end]), A)
    res = typemax(Int)
    for i in 1:length(prex)
        while !isempty(dq) && prex[i] - prex[first(dq)] ≥ K
            res = min(res, i - first(dq))
            popfirst!(dq)
        end
        while !isempty(dq) && prex[i] ≤ prex[last(dq)]
            pop!(dq)
        end
        push!(dq, i)
    end
    res == typemax(Int) ? -1 : res
end
## @lc code=end


text:
	1+1
## --- Setup

    # External packages
    using ProgressMeter: @showprogress
    using Plots
    using JLD, HDF5
    try
        using StatGeochem
    catch
        Pkg.clone("https://github.com/brenhinkeller/StatGeochem.jl")
        using StatGeochem
    end

    # Local utilities
    using HTTP, JSON
    include("RockWeatheringFluxUtilities.jl")

## --- Generate some random points on a sphere

    npoints = 50000;

    rocklat = Array{Float64}(0)
    rocklon = Array{Float64}(0)
    etopo = get_etopo("elevation")
    while length(rocklat) < npoints

        # Generate some random latitudes and longitudes with uniform
        #  spatial density on the globe
        (randlat, randlon) = random_lat_lon(npoints)

        # Find which points are above sea level
        elevations = find_etopoelev(etopo,randlat,randlon)
        abovesea = elevations .> 0

        # Concatenate together all the points that represent exposed crust
        rocklat = vcat(rocklat,randlat[abovesea])
        rocklon = vcat(rocklon,randlon[abovesea])
    end

    rocklat = rocklat[1:npoints]
    rocklon = rocklon[1:npoints]
    elevations = find_etopoelev(etopo,rocklat,rocklon)

## -- Get data from Macrostrat / Burwell API

    zoom = 11
	savefilename = "responses2"
    responses = Array{Any}(length(elevations))
    @showprogress 5 for i = 1:length(elevations)
        try
            responses[i] = query_macrostrat(rocklat[i], rocklon[i], zoom)
        catch
            print("Warning: no data from Macrostrat server \n")
            try
                # Wait and try again
                sleep(10)
                responses[i] = query_macrostrat(rocklat[i], rocklon[i], zoom)
            catch
                # if still nothing, add warnin
                responses[i] = "No response"
                print("Warning: no response from Macrostrat server\n")
            end
        end
        sleep(0.05)
        if mod(i,10000)==0
            save("data/$savefilename.jld", "responses", responses, "elevations", elevations, "latitude", rocklat, "longitude", rocklon, "npoints", npoints)
        end
    end

    save("data/$savefilename.jld", "responses", responses, "elevations", elevations, "latitude", rocklat, "longitude", rocklon, "npoints", npoints)

## --- Load from saved version (if applicable)

    retrive_file = load("data/responses.jld")
    elevations = retrive_file["elevations"]
    rocklat = retrive_file["latitude"]
    rocklon = retrive_file["longitude"]
    responses = retrive_file["responses"]
    npoints = retrive_file["npoints"]

## ---
	retrive_file = load("data/responses2.jld")
	elevations = vcat(elevations, retrive_file["elevations"])
	rocklat = vcat(rocklat, retrive_file["latitude"])
	rocklon = vcat(rocklon, retrive_file["longitude"])
	responses = vcat(responses, retrive_file["responses"])
	npoints = npoints + retrive_file["npoints"]

## --- Parse the macrostrat responses

    # Variables for parsed responses
    rocktype = Array{String}(length(elevations))
    rockdescrip = Array{String}(length(elevations))
    rockname = Array{String}(length(elevations))
    rockstratname = Array{String}(length(elevations))
    rockcomments = Array{String}(length(elevations))
	agemax = Array{Float64}(length(elevations))
	agemin = Array{Float64}(length(elevations))

    # Parse saved responses
    for i = 1:length(elevations)
        rocktype[i] = get_macrostrat_lith(responses[i])
        rockdescrip[i] = get_macrostrat_descrip(responses[i])
        rockname[i] = get_macrostrat_name(responses[i])
        rockstratname[i] = get_macrostrat_strat_name(responses[i])
        rockcomments[i] = get_macrostrat_comments(responses[i])
		agemax[i] = get_macrostrat_max_age(responses[i])
		agemin[i] = get_macrostrat_min_age(responses[i])
    end

    # Convert to lowercase to match names of rock types
    rocktype = lowercase.(rocktype)
    rockdescrip = lowercase.(rockdescrip)
    rockname = lowercase.(rockname)
    rockstratname = lowercase.(rockstratname)
    rockcomments = lowercase.(rockcomments)

    # Replacing tabs with spaces so as not to be confused with the delimitator if exported
    rocktype = replace.(rocktype, "    ", " ")
    rockdescrip = replace.(rockdescrip, "    ", " ")
    rockname = replace.(rockname, "    ", " ")
    rockstratname = replace.(rockstratname, "    ", " ")
    rockcomments = replace.(rockcomments, "    ", " ")

    #= Use this link to check the information on certain points in the macrostrat
    url = "https://macrostrat.org/api/mobile/map_query?lat=$(rocklat[end])&lng=$(rocklon[end])&z=11" =#

## ---

	# Names or partial names of different rock types (from GetBurwellBulkAge.m)

	covertypes = ["cover", "unconsolidated", "quaternary", "lluv", "soil", "regolith", "laterite", "surficial deposits", "talus", "scree", "mass-wasting", "slide", "peat", "swamp", "marsh", "water", "ice", "glaci", "till", "loess", "gravel"] # "debris",
    sedtypes = ["sediment", "fluv", " clast", "siliciclast", "conglomerat", "gravel", "sand", "psamm", "arenit", "arkos", "silt", "mud", "marl", "clay", "shale", "wacke", "argillite", "argillaceous", "pelit", "pebble", "carbonate", "limestone", "dolo", "caliche", "chalk", "travertine", "tavertine", "teravertine", "tufa", "evaporite", " salt", "salt flat", "gypsum", "boulder", "diamict", "tillite", "stream", "beach", "terrace", "chert", "banded iron", "coal", "anthracite", "marine deposits", "turbidite", "flysch", "paleosol"]

    volctypes = ["volcan", "lava", "lahar", "ignimbrite", "ashfall", "tuff", "diatreme", "pipe", "basalt", "andesit", "dacit", "rhyolit", "pillow", "carbonatite", "tephra", "obsidian", "ash", "scoria", "pumice", "cinder", "latite", "basanite", "phonolite", "fonolito", "trachyte", "palagonite", "mugearite", "kimberlite", "ultramafitite", "komatiite",]
	pluttypes = ["pluton", "batholith", "granit", "tonalit", "gabbro", "norite", "diorit", "monzonit", "syenit", "peridot", "dunit", "harzburg", "anorthosite", "mangerite", "charnockite", "pegmatite", "aplite", "trond", "essexite", "pyroxenite", "adamellite", "porphyry", "megacryst", "rapakivi", "bronzitite", "alaskite", "troctolite",]
	hypabyssaltypes = ["intrus", "hypabyssal", "sill", "dike", "stock", "laccolith", "lopolith", "dolerit", "diabase", "porphyry", "microgranite"]

	igntypes = vcat(["igneous", "silicic ", "mafic", "felsic", "basite",],volctypes,pluttypes,hypabyssaltypes)

	metasedtypes = ["para", "metased", "meta-sed", "schist", "quartzite", "marble", "skarn", "slate", "phyllite",]
	metaigntypes = ["ortho", "metaign", "meta-ign", "serpentin", "amphibolit", "greenstone", "eclogite", "metabasite", "migma",]
    mettypes = vcat(metasedtypes, metaigntypes, ["gneiss", "granulit", "hornfels", "granofels", "mylonit", "meta", "cataclasite", "melange", "gouge", "tecton", "calc silicate"])

	lowgradetypes = ["slate", "phyllite", "serpentin", "greenstone", "greenschist", "zeolite", "gossan", "alter", "hydrothermal", "palagonite",]
	highgradetypes = ["crystalline", "basement", "marble", "skarn", "blueschist", "gneiss", "amphibolit", "eclogite", "granulit", "hornfels", "granofels", "sanidinite", "migma", "enderbite", "anorthosite", "charnockite", "pyroxenite", "peridot", "dunit", "harzburg", "high grade metamorphic"]
	cataclastictypes = ["mylonit", "cataclasite", "melange", "gouge", "tecton",]

    # Allocate arrays for each sample for each rock type
    sed = fill(false,npoints)
    ign = fill(false,npoints)
    met = fill(false,npoints)
	volc = fill(false,npoints)
	plut = fill(false,npoints)
	hypabyssal = fill(false,npoints)
	metaign = fill(false,npoints)
	metased = fill(false,npoints)
	lowgrade = fill(false,npoints)
	highgrade = fill(false,npoints)
	cover = fill(false,npoints)

    # Check which burwell "lith" rocktypes match one of the rock types
    # Try the "major:: {...}" type first, if present
    for i = 1:length(sedtypes)
      sed = sed .|  ( match.(r"major.*?{(.*?)}", rocktype) .|> x -> isa(x,RegexMatch) ? contains(x[1], sedtypes[i]) : false )
    end
    for i = 1:length(igntypes)
      ign = ign .|  ( match.(r"major.*?{(.*?)}", rocktype) .|> x -> isa(x,RegexMatch) ? contains(x[1], igntypes[i]) : false )
    end
    for i = 1:length(mettypes)
      met = met .|  ( match.(r"major.*?{(.*?)}", rocktype) .|> x -> isa(x,RegexMatch) ? contains(x[1], mettypes[i]) : false )
    end
	for i = 1:length(covertypes)
	  cover = cover .|  ( match.(r"major.*?{(.*?)}", rocktype) .|> x -> isa(x,RegexMatch) ? contains(x[1], covertypes[i]) : false )
	end

	for i = 1:length(volctypes)
	  volc = volc .|  ( match.(r"major.*?{(.*?)}", rocktype) .|> x -> isa(x,RegexMatch) ? contains(x[1], volctypes[i]) : false )
	end
	for i = 1:length(pluttypes)
      plut = plut .|  ( match.(r"major.*?{(.*?)}", rocktype) .|> x -> isa(x,RegexMatch) ? contains(x[1], pluttypes[i]) : false )
    end
    for i = 1:length(metaigntypes)
      metaign = metaign .|  ( match.(r"major.*?{(.*?)}", rocktype) .|> x -> isa(x,RegexMatch) ? contains(x[1], metaigntypes[i]) : false )
    end
	for i = 1:length(metasedtypes)
	  metased = metased .|  ( match.(r"major.*?{(.*?)}", rocktype) .|> x -> isa(x,RegexMatch) ? contains(x[1], metasedtypes[i]) : false )
	end
	for i = 1:length(lowgradetypes)
      lowgrade = lowgrade .|  ( match.(r"major.*?{(.*?)}", rocktype) .|> x -> isa(x,RegexMatch) ? contains(x[1], lowgradetypes[i]) : false )
    end
	for i = 1:length(highgradetypes)
	  highgrade = highgrade .|  ( match.(r"major.*?{(.*?)}", rocktype) .|> x -> isa(x,RegexMatch) ? contains(x[1], highgradetypes[i]) : false )
	end
	for i = 1:length(hypabyssaltypes)
	  hypabyssal = hypabyssal .|  ( match.(r"major.*?{(.*?)}", rocktype) .|> x -> isa(x,RegexMatch) ? contains(x[1], hypabyssaltypes[i]) : false )
	end

    # Then check the rest of rocktype
    not_matched = .~(sed .| ign .| met .| cover)
    for i = 1:length(sedtypes)
      sed[not_matched] = sed[not_matched] .| contains.(rocktype[not_matched],sedtypes[i])
    end
    for i = 1:length(igntypes)
      ign[not_matched] = ign[not_matched] .| contains.(rocktype[not_matched],igntypes[i])
    end
    for i = 1:length(mettypes)
      met[not_matched] = met[not_matched] .| contains.(rocktype[not_matched],mettypes[i])
    end
	for i = 1:length(covertypes)
      cover[not_matched] = cover[not_matched] .| contains.(rocktype[not_matched],covertypes[i])
    end

	not_matched = .~(sed .| cover .| volc .| plut .| metaign .| metased .| lowgrade .| highgrade)
	for i = 1:length(volctypes)
      volc[not_matched] = volc[not_matched] .| contains.(rocktype[not_matched],volctypes[i])
    end
	for i = 1:length(pluttypes)
      plut[not_matched] = plut[not_matched] .| contains.(rocktype[not_matched],pluttypes[i])
    end
	for i = 1:length(metaigntypes)
	  metaign[not_matched] = metaign[not_matched] .| contains.(rocktype[not_matched],metaigntypes[i])
	end
	for i = 1:length(metasedtypes)
	  metased[not_matched] = metased[not_matched] .| contains.(rocktype[not_matched],metasedtypes[i])
	end
	for i = 1:length(lowgradetypes)
	  lowgrade[not_matched] = lowgrade[not_matched] .| contains.(rocktype[not_matched],lowgradetypes[i])
	end
	for i = 1:length(highgradetypes)
	  highgrade[not_matched] = highgrade[not_matched] .| contains.(rocktype[not_matched],highgradetypes[i])
	end
	for i = 1:length(hypabyssaltypes)
	  hypabyssal[not_matched] = hypabyssal[not_matched] .| contains.(rocktype[not_matched],hypabyssaltypes[i])
	end

    # Then rockname
    not_matched = .~(sed .| ign .| met .| cover)
    for i = 1:length(sedtypes)
      sed[not_matched] = sed[not_matched] .| contains.(rockname[not_matched],sedtypes[i])
    end
    for i = 1:length(igntypes)
      ign[not_matched] = ign[not_matched] .| contains.(rockname[not_matched],igntypes[i])
    end
    for i = 1:length(mettypes)
      met[not_matched] = met[not_matched] .| contains.(rockname[not_matched],mettypes[i])
    end
	for i = 1:length(covertypes)
      cover[not_matched] = cover[not_matched] .| contains.(rockname[not_matched],covertypes[i])
    end

	not_matched = .~(sed .| cover .| volc .| plut .| metaign .| metased .| lowgrade .| highgrade)
	for i = 1:length(volctypes)
      volc[not_matched] = volc[not_matched] .| contains.(rockname[not_matched],volctypes[i])
    end
	for i = 1:length(pluttypes)
      plut[not_matched] = plut[not_matched] .| contains.(rockname[not_matched],pluttypes[i])
    end
	for i = 1:length(metaigntypes)
	  metaign[not_matched] = metaign[not_matched] .| contains.(rockname[not_matched],metaigntypes[i])
	end
	for i = 1:length(metasedtypes)
	  metased[not_matched] = metased[not_matched] .| contains.(rockname[not_matched],metasedtypes[i])
	end
	for i = 1:length(lowgradetypes)
	  lowgrade[not_matched] = lowgrade[not_matched] .| contains.(rockname[not_matched],lowgradetypes[i])
	end
	for i = 1:length(highgradetypes)
	  highgrade[not_matched] = highgrade[not_matched] .| contains.(rockname[not_matched],highgradetypes[i])
	end
	for i = 1:length(hypabyssaltypes)
	  hypabyssal[not_matched] = hypabyssal[not_matched] .| contains.(rockname[not_matched],hypabyssaltypes[i])
	end

    # Then rockdescrip
    not_matched = .~(sed .| ign .| met .| cover)
    for i = 1:length(sedtypes)
      sed[not_matched] = sed[not_matched] .| contains.(rockdescrip[not_matched],sedtypes[i])
    end
    for i = 1:length(igntypes)
      ign[not_matched] = ign[not_matched] .| contains.(rockdescrip[not_matched],igntypes[i])
    end
    for i = 1:length(mettypes)
      met[not_matched] = met[not_matched] .| contains.(rockdescrip[not_matched],mettypes[i])
    end
	for i = 1:length(covertypes)
	  cover[not_matched] = cover[not_matched] .| contains.(rockdescrip[not_matched],covertypes[i])
	end

	not_matched = .~(sed .| cover .| volc .| plut .| metaign .| metased .| lowgrade .| highgrade)
	for i = 1:length(volctypes)
	  volc[not_matched] = volc[not_matched] .| contains.(rockdescrip[not_matched],volctypes[i])
	end
	for i = 1:length(pluttypes)
	  plut[not_matched] = plut[not_matched] .| contains.(rockdescrip[not_matched],pluttypes[i])
	end
	for i = 1:length(metaigntypes)
	  metaign[not_matched] = metaign[not_matched] .| contains.(rockdescrip[not_matched],metaigntypes[i])
	end
	for i = 1:length(metasedtypes)
	  metased[not_matched] = metased[not_matched] .| contains.(rockdescrip[not_matched],metasedtypes[i])
	end
	for i = 1:length(lowgradetypes)
	  lowgrade[not_matched] = lowgrade[not_matched] .| contains.(rockdescrip[not_matched],lowgradetypes[i])
	end
	for i = 1:length(highgradetypes)
	  highgrade[not_matched] = highgrade[not_matched] .| contains.(rockdescrip[not_matched],highgradetypes[i])
	end
	for i = 1:length(hypabyssaltypes)
	  hypabyssal[not_matched] = hypabyssal[not_matched] .| contains.(rockdescrip[not_matched],hypabyssaltypes[i])
	end

	# Exclude suspected cover from the other three categories, just to be sure
	ign = ign .& .~cover
	met = met .& .~cover
	sed = sed .& .~cover
	cryst = ign .| (met .& .~sed)

    not_matched = .~(sed .| ign .| met .| cover);
    multi_matched = (sed .& ign) .| (sed .& met) .| (ign .& met);

    number_not_matched = sum(not_matched)
    number_multi_matched = sum(multi_matched)

    print("not matched = $number_not_matched, conflicting matches = $number_multi_matched\n")

## ---

    # Print proportions of sed vs ign, met, and nonsed rocks in responses
    print("sed = ", sum(sed), ", ign + (met & ~sed)= ", sum(ign .| (met .& .~sed)))

    # Create a file to check matching errors
    writedlm("notmatched.tsv", hcat(rocktype[not_matched], rockname[not_matched], rockdescrip[not_matched], rockstratname[not_matched], rockcomments[not_matched]))

    writedlm("multimatched.tsv", hcat(rocktype[multi_matched], rockname[multi_matched], rockdescrip[multi_matched], rockstratname[multi_matched], rockcomments[multi_matched]))

    writedlm("ignsed.tsv", hcat(rocktype[ign .& sed], rockname[ign .& sed], rockdescrip[ign .& sed], rockstratname[ign .& sed], rockcomments[ign .& sed]))

	writedlm("plutonic.tsv", hcat(rocktype[plut], rockname[plut], rockdescrip[plut], rockstratname[plut], rockcomments[plut]))

	writedlm("intrusive.tsv", hcat(rocktype[hypabyssal .& .~plut], rockname[hypabyssal .& .~plut], rockdescrip[hypabyssal .& .~plut], rockstratname[hypabyssal .& .~plut], rockcomments[hypabyssal .& .~plut]))

	t = elevations .> 4000
	writedlm("highelev.tsv", hcat(rocktype[t], rockname[t], rockdescrip[t], rockstratname[t], rockcomments[t]))

## --- Plot abundance of sed, highgrade+plutonic and lowgrade+volcanic

	nresampled = 10^6
	edges = 0:40:4000
	centers = cntr(edges)
	nbins = length(centers)

	t = sed .& .~cover
	age = (agemax[t] + agemin[t])/2
	# sigma = (agemax[t] - agemin[t])/2/sqrt(2*log(2))
	sigma = (agemax[t] - agemin[t])/2
	# ageresampled = bsresample_unif(age,sigma,nresampled)
	ageresampled = bsresample_unif_norm(age,sigma,age*0.025,nresampled)
	ageresampled = ageresampled[(ageresampled .> 0) .& (ageresampled .< 4500)]
	h = histogram(ageresampled, bins=edges, fill=(0,:darkblue), linewidth=0.1, label="")
	plot!(h,ylabel="N", ylims=(0,nresampled/nbins*6), xlims=(0,4000), xflip=true, xlabel="Age (Ma)")
	savefig(h,"sed.pdf")

	t = ign .& .~(cover .| sed)
	age = (agemax[t] + agemin[t])/2
	sigma = (agemax[t] - agemin[t])/2/sqrt(2*log(2))
	# ageresampled = bsresample_unif(age,sigma,nresampled)
	ageresampled = bsresample_unif_norm(age,sigma,age*0.025,nresampled)
	ageresampled = ageresampled[(ageresampled .> 0) .& (ageresampled .< 4500)]
	h = histogram(ageresampled, bins=edges, fill=(0,:darkred), linewidth=0.1, label="")
	plot!(h, ylabel="N", ylims=(0,nresampled/nbins*6), xlims=(0,4000), xflip=true,  xlabel="Age (Ma)")
	savefig(h,"ign.pdf")

	t = (volc .| lowgrade) .& .~(cover .| sed)
	age = (agemax[t] + agemin[t])/2
	sigma = (agemax[t] - agemin[t])/2/sqrt(2*log(2))
	ageresampled = bsresample_unif(age,sigma,nresampled)
	ageresampled = ageresampled[(ageresampled .> 0) .& (ageresampled .< 4500)]
	h = histogram(ageresampled, bins=edges, fill=(0,:red), linewidth=0.1, label="")
	plot!(h, ylabel="N", ylims=(0,nresampled/nbins*6), xlims=(0,4000), xflip=true,  xlabel="Age (Ma)")
	savefig(h,"volc.pdf")

	t = hypabyssal .& .~(cover .| sed)
	age = (agemax[t] + agemin[t])/2
	sigma = (agemax[t] - agemin[t])/2/sqrt(2*log(2))
	# ageresampled = bsresample_unif(age,sigma,nresampled)
	ageresampled = bsresample_unif_norm(age,sigma,age*0.025,nresampled)
	ageresampled = ageresampled[(ageresampled .> 0) .& (ageresampled .< 4500)]
	h = histogram(ageresampled, bins=edges, fill=(0,:darkred), linewidth=0.1, label="")
	plot!(h, ylabel="N", ylims=(0,nresampled/nbins*6), xlims=(0,4000), xflip=true,  xlabel="Age (Ma)")
	savefig(h,"shallow intrusive.pdf")

	t = plut .& .~(cover .| sed)
	age = (agemax[t] + agemin[t])/2
	sigma = (agemax[t] - agemin[t])/2/sqrt(2*log(2))
	# ageresampled = bsresample_unif(age,sigma,nresampled)
	ageresampled = bsresample_unif_norm(age,sigma,age*0.025,nresampled)
	ageresampled = ageresampled[(ageresampled .> 0) .& (ageresampled .< 4500)]
	h = histogram(ageresampled, bins=edges, fill=(0,:darkred), linewidth=0.1, label="")
	plot!(h, ylabel="N", ylims=(0,nresampled/nbins*6), xlims=(0,4000), xflip=true,  xlabel="Age (Ma)")
	savefig(h,"plut.pdf")

	t = highgrade .& .~(cover .| sed)
	age = (agemax[t] + agemin[t])/2
	sigma = (agemax[t] - agemin[t])/2/sqrt(2*log(2))
	# ageresampled = bsresample_unif(age,sigma,nresampled)
	ageresampled = bsresample_unif_norm(age,sigma,age*0.025,nresampled)
	ageresampled = ageresampled[(ageresampled .> 0) .& (ageresampled .< 4500)]
	h = histogram(ageresampled, bins=edges, fill=(0,:darkred), linewidth=0.1, label="")
	plot!(h, ylabel="N", ylims=(0,nresampled/nbins*6), xlims=(0,4000), xflip=true,  xlabel="Age (Ma)")
	savefig(h,"highgrade.pdf")

	t = (plut .| highgrade) .& .~(cover .| sed)
	age = (agemax[t] + agemin[t])/2
	sigma = (agemax[t] - agemin[t])/2/sqrt(2*log(2))
	# ageresampled = bsresample_unif(age,sigma,nresampled)
	ageresampled = bsresample_unif_norm(age,sigma,age*0.025,nresampled)
	ageresampled = ageresampled[(ageresampled .> 0) .& (ageresampled .< 4500)]
	h = histogram(ageresampled, bins=edges, fill=(0,:darkred), linewidth=0.1, label="")
	plot!(h, ylabel="N", ylims=(0,nresampled/nbins*6), xlims=(0,4000), xflip=true,  xlabel="Age (Ma)")
	savefig(h,"basement (plut+highgrade).pdf")

	t = .~(cover .| sed)
	age = (agemax[t] + agemin[t])/2
	sigma = (agemax[t] - agemin[t])/2/sqrt(2*log(2))
	# ageresampled = bsresample_unif(age,sigma,nresampled)
	ageresampled = bsresample_unif_norm(age,sigma,age*0.025,nresampled)
	ageresampled = ageresampled[(ageresampled .> 0) .& (ageresampled .< 4500)]
	h = histogram(ageresampled, bins=edges, fill=(0,:darkred), linewidth=0.1, label="")
	plot!(h, ylabel="N", ylims=(0,nresampled/nbins*6), xlims=(0,4000), xflip=true,  xlabel="Age (Ma)")
	savefig(h,"nonsed.pdf")

	t = .~(cover .| sed .| metased .| volc .| hypabyssal .| not_matched)
	age = (agemax[t] + agemin[t])/2
	sigma = (agemax[t] - agemin[t])/2/sqrt(2*log(2))
	# ageresampled = bsresample_unif(age,sigma,nresampled)
	ageresampled = bsresample_unif_norm(age,sigma,age*0.025,nresampled)
	ageresampled = ageresampled[(ageresampled .> 0) .& (ageresampled .< 4500)]
	h = histogram(ageresampled, bins=edges, fill=(0,:darkred), linewidth=0.1, label="")
	plot!(h, ylabel="N", ylims=(0,nresampled/nbins*6), xlims=(0,4000), xflip=true,  xlabel="Age (Ma)")
	savefig(h,"nonsednonvolc.pdf")

## --- get references

	authors = get_macrostrat_ref_authors.(responses)
	years = get_macrostrat_ref_year.(responses)
	titles = get_macrostrat_ref_title.(responses)
	dois = get_macrostrat_ref_doi.(responses)

	refstrings = authors .* "|" .* years .* "|" .* titles .* "|" .* dois
	writedlm("burwellrefs.tsv",unique(refstrings),"\t")

## --- Find slope and erosion rate

    function Emmkyr(slp)
        return 10^(slp*0.00567517 + 0.971075)
    end

    slope = get_srtm15plus_aveslope("slope")
    rockslope = find_srtm15plus_aveslope_around(slope,rocklat,rocklon, halfwidth=7, max_allowed_slope=0xffff)

	# Plot slope-area histogram
	h = histogram(rockslope, label="",xlims=[0,500])
	plot!(h,xlabel="Slope (m/km)",ylabel="N")
	savefig(h,"Slope_Histogram.pdf")
	display(h)

    rockEmmkyr = Emmkyr.(rockslope)

    sedEsum = nansum(rockEmmkyr[sed])
    crystEsum = nansum(rockEmmkyr[cryst])

    sedEmean = nanmean(rockEmmkyr[sed])
    crystEmean = nanmean(rockEmmkyr[cryst])
	Emean = nanmean(rockEmmkyr)

    propotionSedArea = sum(sed)/(sum(sed)+sum(ign .| (met .& .~sed)));
    print("$(round(100*propotionSedArea,3))% sed area, $(round(100*(1-propotionSedArea),3))% crystalline area\n")
    print("sed E sum: $sedEsum cryst E sum: $crystEsum\n")
    print("$(round(100*sedEsum/(sedEsum+crystEsum),3))% sed flux, $(round(100*crystEsum/(sedEsum+crystEsum),3))% crystalline flux\n")
    print("sed mean E: $(round(sedEmean,3)), crystalline mean E: $(round(crystEmean,3)), overall mean: $(round(Emean,3))\n")

## --- Proportion vs slope

nbins = 50
xmin = 0
xmax = 500

prop_cryst = Array{Float64}(nbins)
prop_sed = Array{Float64}(nbins)
binedges = range(xmin,xmax,length=nbins+1)
bincenters = cntr(binedges)
for i=1:nbins
	t = (rockslope.>binedges[i]) .& (rockslope.<binedges[i+1])
	prop_cryst[i] = sum(cryst[t])/sum(cryst[t] .| sed[t])
	prop_sed[i] = sum(sed[t])/sum(cryst[t] .| sed[t])
end

h = plot(bincenters,prop_sed,label = "", xlabel="Slope (m/km)", ylabel="Proportion of Sedimentary Bedrock")
# plot!(h,bincenters,prop_cryst+prop_sed,label="Sedimentary Basement")
plot!(legend=:topleft)
savefig(h,"Proportion_sedimentary_bedrock_slope.pdf")
display(h)

## --- Proportion vs elevation

nbins = 50
xmin = 0
xmax = 5000

prop_cryst = Array{Float64}(nbins)
prop_sed = Array{Float64}(nbins)
binedges = range(xmin,xmax,length=nbins+1)
bincenters = cntr(binedges)
for i=1:nbins
	t = (elevations.>binedges[i]) .& (elevations.<binedges[i+1])
	prop_cryst[i] = sum(cryst[t])/sum(cryst[t] .| sed[t])
	prop_sed[i] = sum(sed[t])/sum(cryst[t] .| sed[t])
end

h = plot(bincenters,prop_sed,label = "", xlabel="Elevation (m)", ylabel="Proportion of Sedimentary Bedrock")
# plot!(h,bincenters,prop_cryst+prop_sed,label="Sedimentary Basement")
savefig(h,"Proportion_sedimentary_bedrock_elevation.pdf")
display(h)
## --- For comparison

    # Phanerozoic
    sedimentSubduction_km3_yr = 1.64 # Clift, 2009
    continentSedimentation_km3_yr = 0.8616 # Macrostrat
    continentalErosion_km3_yr = sedimentSubduction_km3_yr + sedimentSubduction_km3_yr

    continentalErosionQuaternary_km3_yr = 6.5 #\cite{Clift:2009fm}  / Milliman (1990)
    continentalArea_km2 =  1.4894E8 # km^2

    continentalErosion_km_yr = continentalErosion_km3_yr / continentalArea_km2
    continentalErosion_mm_kyr = continentalErosion_km_yr * 1000 * 1000 * 1000

### ---
    using StatGeochem
    elev = get_srtm15plus("elevation")

    function movmean33(x::Array{<:Any,2})
        # Simple moving average of 9-element box for 2-d array
        m = Array{Float64}(size(x))
        for i = 2:size(x,1)-1
            for j = 2:size(x,2)-1
                adj = [x[i,j],x[i+1,j+1],x[i,j+1],x[i-1,j+1],x[i-1,j],x[i-1,j-1],x[i,j-1],x[i+1,j-1],x[i+1,j]]
                m[i,j] = mean(adj)
            end
        end
        m[1,:] = m[2,:]
        m[end,:] = m[end-1,:]
        m[:,1] = m[:,2]
        m[:,end] = m[:,end-1]
        return m;
    end

    elev_smoothed = movmean33(elev)

    # Save results
    fid = h5open("srtm15plus_elev_smoothed.h5","w");
    g = g_create(fid, "vars");
    print("Saving to HDF5:\n")
    g["y_lat_cntr"] = y_lat_cntr
    g["x_lon_cntr"] = x_lon_cntr
    g["cellsize"] = cellsize
    g["scalefactor"] = scalefactor
    @time g["elevation","compress",3] = elev_smoothed;
    close(fid)

### ---
	elev = elev_smoothed
	inbasin = fill(false,size(elev))

	inbasin[elev .< -100] = true

	function fill_uphill(elev,inbasin,i,j,depth=0,tolerance=15)

        inbasin[i,j] = true
        depth += 1

        # If we're about to run in the recursion limit or the edge of the map, stop
        if depth > 10000 || i<2 || i > size(elev,1)-1 || j<2 || j > size(elev,2)-1
            return inbasin
        end

		# Otherwise, check each adjoining square recursively
		e = elev[i,j] - tolerance

		if elev[i+1,j+1] >= e && ~inbasin[i+1,j+1]
			inbasin = fill_uphill(elev,inbasin,i+1,j+1,depth)
		end
		if elev[i,j+1] >= e && ~inbasin[i,j+1]
			inbasin = fill_uphill(elev,inbasin,i,j+1,depth)
		end
		if elev[i-1,j+1] >= e && ~inbasin[i-1,j+1]
			inbasin = fill_uphill(elev,inbasin,i-1,j+1,depth)
		end
		if elev[i-1,j] >= e && ~inbasin[i-1,j]
			inbasin = fill_uphill(elev,inbasin,i-1,j,depth)
		end
		if elev[i-1,j-1] >= e && ~inbasin[i-1,j-1]
			inbasin = fill_uphill(elev,inbasin,i-1,j-1,depth)
		end
		if elev[i,j-1] >= e && ~inbasin[i,j-1]
			inbasin = fill_uphill(elev,inbasin,i,j-1,depth)
		end
		if elev[i+1,j-1] >= e && ~inbasin[i+1,j-1]
			inbasin = fill_uphill(elev,inbasin,i+1,j-1,depth)
		end
        if elev[i+1,j] >= e && ~inbasin[i+1,j]
            inbasin = fill_uphill(elev,inbasin,i+1,j,depth)
        end
		return inbasin
	end

	@showprogress for i=2:size(inbasin,1)-1
		for j=2:size(inbasin,2)-1
			if inbasin[i,j]
                # inbasin = fill_uphill(elev,inbasin,i,j)
				if ~inbasin[i+1,j] || ~inbasin[i+1,j+1] || ~inbasin[i,j+1] || ~inbasin[i-1,j+1] || ~inbasin[i-1,j] || ~inbasin[i-1,j-1]  || ~inbasin[i,j-1]  || ~inbasin[i+1,j-1]
					inbasin = fill_uphill(elev,inbasin,i,j)
				end
			end
		end
	end

    heatmap(downsample(inbasin,10))

    1+1


text:
# MIT License
#
# Copyright (c) 2018 Martin Biel
#
# Permission is hereby granted, free of charge, to any person obtaining a copy
# of this software and associated documentation files (the "Software"), to deal
# in the Software without restriction, including without limitation the rights
# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
# copies of the Software, and to permit persons to whom the Software is
# furnished to do so, subject to the following conditions:
#
# The above copyright notice and this permission notice shall be included in all
# copies or substantial portions of the Software.
#
# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
# SOFTWARE.

@with_kw mutable struct BBParameters{T <: AbstractFloat}
    γ₀::T = 0.1
end

"""
    BBStep

Functor object for using the Barzilai-Borwein step size in a quasigradient algorithm. Create by supplying a [`BB`](@ref) object through `step` to `QuasiGradient.Optimizer` or by setting the [`StepSize`](@ref) attribute.

...
# Parameters
- `γ₀::AbstractFloat = 0.1`: Initial step-size and fallback if BB is numerically unstable
...
"""
struct BBStep{T <: AbstractFloat} <: AbstractStep
    parameters::BBParameters{T}
    xprev::Vector{T}
    gprev::Vector{T}

    function BBStep(::Type{T}; kw...) where T <: AbstractFloat
        return new{T}(BBParameters{T}(; kw...), Vector{T}(), Vector{T}())
    end
end

function step(step::BBStep, k::Integer, fval::Float64, x::AbstractVector, g::AbstractVector)
    if k == 1
        resize!(step.xprev, length(x))
        resize!(step.gprev, length(x))
        step.xprev .= x
        step.gprev .= g
        return step.parameters.γ₀
    end
    s = x - step.xprev
    y = g - step.gprev
    sy = (s⋅y)
    if abs(sy) <= sqrt(eps())
        return step.parameters.γ₀
    end
    η = (norm(s,2) ^ 2) / sy
    step.xprev .= x
    step.gprev .= g
    return η
end

# API
# ------------------------------------------------------------
"""
    BB

Factory object for [`BBStep`](@ref). Pass to `step` in `Quasigradient.Optimizer` or set the [`StepSize`](@ref) attribute. See ?BBStep for parameter descriptions.

"""
struct BB <: AbstractStepSize
    parameters::BBParameters{Float64}
end
BB(γ₀::AbstractFloat) = BB(BBParameters(; γ₀ = Float64(γ₀)))
BB(; kw...) = BB(BBParameters(; kw...))

function (stepsize::BB)(::Type{T}) where T <: AbstractFloat
    return BBStep(T; type2dict(stepsize.parameters)...)
end

function str(::BB)
    return "BB step size"
end


text:
"""
    get_results(model, case)

Extract the results of a simulation, store in a struct.
"""
function get_results(f::Float64, case::Case, demand_flexibility::DemandFlexibility)::Results
    status = "OPTIMAL"
    sets = _make_sets(case; storage=nothing, demand_flexibility=demand_flexibility)
    # These variables will always be in the results
    pg = JuMP.value.(m[:pg])
    pf = JuMP.value.(m[:pf])
    lmp = -1 * JuMP.shadow_price.(m[:powerbalance])
    congl_temp = -1 * JuMP.shadow_price.(m[:branch_min])
    congu_temp = -1 * JuMP.shadow_price.(m[:branch_max])
    # If DC lines are present, separate their results
    # Initialize with empty arrays, to be discarded later if they stay empty
    pf_dcline = zeros(0, 0)
    num_dclines = length(case.dclineid)
    if num_dclines > 0
        pf_dcline = pf[(end - num_dclines + 1):end, :]
        pf = pf[1:(end - num_dclines), :]
    end
    # Ensure that we report congestion on all branches, even infinite capacity
    num_hour = size(pf, 2)
    congl = zeros(sets.num_branch_ac, num_hour)
    congu = zeros(sets.num_branch_ac, num_hour)
    # Access congl_temp via key `i`, then store result in congl at position `i`
    for i in intersect(Set(sets.noninf_branch_idx), Set(1:(sets.num_branch_ac)))
        congl[i, :] = congl_temp[i, :]
        congu[i, :] = congu_temp[i, :]
    end
    # These variables will only be in the results if the model has storage
    # Initialize with empty arrays, to be discarded later if they stay empty
    storage_pg = zeros(0, 0)
    storage_e = zeros(0, 0)
    try
        storage_dis = JuMP.value.(m[:storage_dis])
        storage_chg = JuMP.value.(m[:storage_chg])
        storage_e = JuMP.value.(m[:storage_soc])
        storage_pg = storage_dis - storage_chg
    catch e
        if isa(e, KeyError)
            # Thrown when storage variables are not defined in the model
        else
            # Unknown error, rethrow it
            rethrow(e)
        end
    end

    # This variable will only be in the results if load shedding is enabled
    # Initialize with empty arrays, to be discarded later if they stay empty
    load_shed = zeros(0, 0)
    try
        load_shed_temp = JuMP.value.(m[:load_shed])
        load_shed = sets.load_bus_map * load_shed_temp
    catch e
        if isa(e, KeyError)
            # Thrown when load_shed is not defined in the model
        else
            # Unknown error, rethrow it
            rethrow(e)
        end
    end

    # These variables will only be in the results if the model has flexible demand
    # Initialize with empty arrays, to be discarded later if they stay empty
    load_shift_up = zeros(0, 0)
    load_shift_dn = zeros(0, 0)
    try
        load_shift_up_temp = JuMP.value.(m[:load_shift_up])
        load_shift_dn_temp = JuMP.value.(m[:load_shift_dn])
        load_shift_up = sets.flexible_load_bus_map * load_shift_up_temp
        load_shift_dn = sets.flexible_load_bus_map * load_shift_dn_temp
    catch e
        if isa(e, KeyError)
            # Thrown when load shift variables are `nothing`
        else
            # Unknown error, rethrow it
            rethrow(e)
        end
    end

    trans_viol = zeros(0, 0)
    try
        trans_viol = JuMP.value.(m[:trans_viol])
    catch e
        if isa(e, MethodError)
            # Thrown when trans_viol is `nothing`
        elseif isa(e, KeyError)
            # Thrown when trans_viol is not defined in the model
        else
            # Unknown error, rethrow it
            rethrow(e)
        end
    end

    results = Results(;
        pg=pg,
        pf=pf,
        lmp=lmp,
        congl=congl,
        congu=congu,
        pf_dcline=pf_dcline,
        f=f,
        storage_pg=storage_pg,
        storage_e=storage_e,
        load_shed=load_shed,
        load_shift_up=load_shift_up,
        load_shift_dn=load_shift_dn,
        status=status,
        trans_viol=trans_viol,
    )
    return results
end


text:
using AoGExtensions
using Documenter

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

makedocs(;
    modules=[AoGExtensions],
    authors="Fabian Greimel <fabgrei@gmail.com> and contributors",
    repo="https://github.com/greimel/AoGExtensions.jl/blob/{commit}{path}#{line}",
    sitename="AoGExtensions.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://greimel.github.io/AoGExtensions.jl",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
    ],
)

deploydocs(;
    repo="github.com/greimel/AoGExtensions.jl",
)


text:
export clustertree

#
# clustertree partitions a collection x by maximizing distances normalized by the
# minimal diameter along coordinate axes in 2D. It uses complex numbers to represent
# points in the plane R². The Admissibility condition for low-rank approximants to
# asymptotically smooth kernels is:
#
# Adm_η(σ×τ) = min{diam(Γ_σ),diam(Γ_τ)} ≤ 2η dist(Γ_σ,Γ_τ),
#
# and clustertree creates partitions that minimize η, the strength of the admissibility.
#

clustertree(x) = x
clustertree(sp::PiecewiseSpace) = Space(clustertree(domain(sp)))
clustertree(d::UnionDomain) = clustertree(d.domains)

clustertree(x::Tuple) = clustertree(collect(x))

function clustertree(x::AbstractVector)
    A = weightedadjacency(x)
    c = map(centroid,x)
    clustertree(A,c,x)
end

clustertree(A::AbstractMatrix,c::AbstractVector,x::AbstractVector) = clustertree(A,c,x,1:size(A,1))

for (TYP,HTYP) in ((:Number,:HierarchicalVector),(:Domain,:HierarchicalDomain))
    @eval begin
        function clustertree{T<:$TYP}(A::AbstractMatrix,c::AbstractVector,x::AbstractVector{T},p::AbstractVector)
            i,j = indmaxp(A,p)

            cmax = c[i]-c[j]
            rbar,ibar = abs2(real(cmax)),abs2(imag(cmax))
            rmid,imid = reim((c[i]+c[j])/2)

            (f,mid) = rbar ≥ ibar ? (real,rmid) : (imag,imid)

            np₁,np₂ = 0,0
            for k in p
                f(c[k]) < mid ? np₁+=1 : np₂+=1
            end
            p₁,p₂ = zeros(Int,np₁),zeros(Int,np₂)
            kp₁,kp₂ = 1,1
            for k in p
                f(c[k]) < mid ? (p₁[kp₁] = k; kp₁+=1) : (p₂[kp₂] = k; kp₂+=1)
            end
            np₁ == 0 && (push!(p₁,p₂[1]); shift!(p₂); np₁+=1; np₂-=1)
            np₂ == 0 && (push!(p₂,p₁[1]); shift!(p₁); np₁-=1; np₂+=1)

            if np₁ == 1 && np₂ == 1
                return $HTYP((x[p₁[1]],x[p₂[1]]))
            elseif np₁ == 1
                return $HTYP((x[p₁[1]],clustertree(A,c,x,p₂)))
            elseif np₂ == 1
                return $HTYP((clustertree(A,c,x,p₁),x[p₂[1]]))
            else
                return $HTYP((clustertree(A,c,x,p₁),clustertree(A,c,x,p₂)))
            end
        end
    end
end

#
# The weighted adjacency matrix computes all the possible distances normalized by the
# minimum of the diameters of the sets. This is related to the Admissibility condition.
#

function weightedadjacency(x::AbstractVector)
    n = length(x)
    A = zeros(real(mapreduce(eltype,promote_type,x)),n,n)
    for j=1:n,i=1:j
        @inbounds A[i,j] = dist2(x[i],x[j])/min(diam2(x[i]),diam2(x[j]))
    end
    for j=1:n,i=1:j-1
        @inbounds A[j,i] = A[i,j]
    end
    A
end

# find the maximum of A over the indices p without allocating A[p,p].

function findmaxp(a,p)
    if isempty(a)
        throw(ArgumentError("collection must be non-empty"))
    end
    m = a[p[1],p[1]]
    mij = p[1],p[1]
    for i in p,j in p
        aij = a[i,j]
        if aij > m || m!=m
            m = aij
            mij = i,j
        end
    end
    return (m, mij)
end
indmaxp(a,p) = findmaxp(a,p)[2]


text:
# to run tests, type in Julia: Pkg.test("ModalParam")

include("../src/modal_test.jl")


text:
# Core data structures for analyzing Wordle

# const names that require reading files or expensive computations are not recomputed.

import Base: empty!, getindex, length, show

"""
    Return vector of strings from a file.
    Canonicalize uppercase as lowercase.
"""
readWords(filename::AbstractString) = String[lowercase(line) for line in eachline(filename)]

if !isdefined(Main, :answerWords)
"""
    Possible answer words
"""
const answerWords = sort(readWords("answers.txt"))
end

"""
     Index into answerWords
"""
const AnswerIndex = UInt16

if !isdefined(Main, :guessWords)
"""
    Allowed guess words
"""
const guessWords = sort([answerWords; readWords("guesses-other.txt")])
end

"""
     Index into guessWords
"""
const GuessIndex = UInt16

"""
     Get index of word in list

Example:
```
julia> indexOfWord("abate",answerWords)
0x0003
```
"""
function indexOfWord(word::String, list::Vector{String})
    i = searchsortedfirst(list, word)
    if i > length(list) || list[i] ≠ word
        error("$word not in list")
    end
    return UInt16(i)
end

"""
     Get indices of word in list
"""
indicesOfWords(words::Vector{String}, list::Vector{String}) = [indexOfWord(w, list) for w in words]

"""
    Color code for a position in a response
"""
@enum Color black=0 yellow=1 green=2

"""
    Color enums for letters
"""
const colorOfLetter = Dict('b' => black, 'y' => yellow, 'g' => green)

"""
    Response to a guess.

    Can be constructed from string of letter codes or (answer, guess).

Examples:
```
julia> Response("bbbyg")
BBBYG

julia> Response("abbot","abate")
GGYYB
```
"""
struct Response
    rep :: UInt8 # Encoded in ternary, with 1 added so it can be used as a one-based index.

    function Response(ternaryEncoding::Integer)
        @assert ternaryEncoding ≥ 1
        @assert ternaryEncoding ≤ 3^5
        return new(ternaryEncoding)
    end

    function Response(colors::AbstractString)
        @assert length(colors) == 5
        sum = 1
        placeValue = 1
        for c in colors
            sum += Int(colorOfLetter[lowercase(c)]) * placeValue
            placeValue *= 3
        end
        return new(sum)
    end

    function Response(answer::String, guess::String)
        @assert length(answer) == 5
        @assert length(guess) == 5
        sum = 1
        placeValue = 1
        for i in 1:5
            if guess[i] == answer[i]
                sum += Int(green) * placeValue
            elseif guess[i] in answer
                sum += Int(yellow) * placeValue
            end
            placeValue *= 3
        end
        return new(sum)
    end
end

"""
    Print response using colored letters.
"""
function show(io::IO, response::Response)
    value = response.rep - 1
    style = [:light_black, :light_yellow, :light_green]
    for position in 1:5
        k = value % 3
        printstyled(io, "byg"[k+1]; color=style[k+1])
        value = div(value, 3)
    end
end

if !isdefined(Main, :responseMatrix)
"""
    Matrix of precomputed responses, indexed by [answerIndex, guessIndex].
"""
const responseMatrix = [Response(answerWords[i], guessWords[j]) for i in 1:length(answerWords), j in 1:length(guessWords)]
end

"""
    Partition

Partition of answer indices. The interface behaves similarly to an assocative map from Response to Vector{AnswerIndex}.
"""
struct Partition
    # Indices of answers in the bucket, or nothing if bucket not yet needed.
    buckets::Vector{Union{Nothing,Vector{AnswerIndex}}}

    # Indices of non-empty buckets
    nonEmptyBuckets::Vector{Int16}

    # Empty partition
    function Partition()
        emptyBuckets = [Int16[] for i in 1:3^5]
        return new(fill(nothing, 3^5), Int16[])
    end
end

"""
    Support for iterating over non-empty buckets in a partition. Each iteration returns a (Response, Vector{AnswerIndex}).
"""
function Base.iterate(p::Partition, state=1)
    if length(p.nonEmptyBuckets) < state
        return nothing
    else
        i = p.nonEmptyBuckets[state]
        return ((Response(i), p.buckets[i]), state + 1)
    end
end

"""
    Number of non-empty buckets in a partition
"""
length(p::Partition) = length(p.nonEmptyBuckets)

"""
    Get bucket for given response
"""
function getindex(p::Partition, r::Response)
    b = p.buckets[r.rep]
    if b == nothing
        b = AnswerIndex[]
        p.buckets[r.rep] = b
    end
    return b
end

function show(io::IO, p::Partition)
    rowDelim = ""
    for (r, b) in p
        print(io, rowDelim)
        rowDelim = "\n"
        print(io, r)
        for j in b
            print(io, " ", answerWords[j])
        end
    end
end

"""
    Set all buckets to empty.
"""
function empty!(p::Partition)
    for i in p.nonEmptyBuckets
        empty!(p.buckets[i])
    end
    empty!(p.nonEmptyBuckets)
end

"""
    Add answer to partition
"""
function pushAt!(p::Partition, r::Response, a::AnswerIndex)
    b = p[r]
    if isempty(b)
        push!(p.nonEmptyBuckets, r.rep)
    end
    push!(b, a)
    return p
end

"""
    Set p to partition of answers induced by guess.
    This is a low-level imperative method that takes indices instead of strings.
"""
function partition!(p::Partition, answers::Vector{AnswerIndex}, guess::GuessIndex)
    empty!(p)
    for a in answers
        pushAt!(p, responseMatrix[a, guess], a)
    end
end

"""
    partition(answers, guess)

Return partition of answers induced by guess.
This is a high-level method that takes strings instead of indices.

Example:
```
julia> partition(["aback", "abase", "abate", "abbey"], "abbot")
ggybb aback abase
ggyby abate
gggbb abbey
```
"""
function partition(answers::Vector{String}, guess::String)
    p = Partition()
    partition!(p, indicesOfWords(answers, answerWords), indexOfWord(guess, guessWords))
    return p
end

"""
    partition[response]

Returns answers in partition for given response.
This is a high-level method that takes the string form of a response and returns strings.
```
julia> partition(["aback", "abase", "abate", "abbey"], "abbot")["ggybb"]
2-element Vector{String}:
 "aback"
 "abase"
```
"""
getindex(p::Partition, response::String) = answerWords[p[Response(response)]]

"""
    winnow(answers, guess, response)

Return subset of answers consistent with response to a guess.
"""
function winnow(answers::Vector{String}, guess::String, response::String)
    r = Response(response)
    g = indexOfWord(guess, guessWords)
    a = indicesOfWords(answers, answerWords)
    return answerWords[[ai for ai in a if responseMatrix[ai, g] == r]]
end

nothing

max_stars_repo_name: mbeltagy/WordleAnalysis

max_issues_repo_name: mbeltagy/WordleAnalysis

max_forks_repo_name: mbeltagy/WordleAnalysis


text:
include("../particlePlot/particlePlot.jl")

# State x = [pₓ, vₓ , ϕ₁, ϕ₂, β₁, β₂, r₁, r₂, a]
# Noise u = [a, ϕ₁, ϕ₂, β₁, β₂]
# Out   y₁ = ∑rₖe^(iβₖ)e^(i2π<pₖ,u(ϕₖ)>)+ny
#       y₂ = a + na
#
# function runTest(N, method, debug, T = 20, withImprovement=true)
#   if withImprovement
#     f! = signalUpdate!
#   else
#     f! = signalUpdateOld!
#   end
#   #g = signalOutNoNoise
#   g = signalOutNoNoiseFast
#   gn = signalOut
#   Ny = 3
#   ginv = signalpxy
#   randsf = (randGen!, :aux)
#   #est = (x,w,t) -> x[:,findmax(w)[2]]
#   est = (x,w,t) -> sum(x.*exp(w'),2)
#
#   x0 = [0., 0, .2, .7, 0., 0. , .8, 1.2, 0]
#   x, y = generateRealSequence(f!, gn, x0, T, randn!, Ny, Nu = 8)
#   xhat = method(f!, g, ginv, randsf, x0, y, est, N = N, Nu = 5, debug=debug, xreal=x)
#   x, xhat
# end
#
function rms(x)
  sqrt(1/length(x)*sum(x.^2))
end

function testRadio()
  debug = true
  Ns = 2.^(5:5)
  M = 1
  T = 100
  RMS = Array{Float64}(length(Ns),M,2)
  largeRMS = Array{Float64}(length(Ns),2)
  rmsMean = Array{Float64,2}(length(Ns),2)
  rmsVariance = Array{Float64,2}(length(Ns),2)
  # Method, Particles, Run nbr, variable, time
  allxhats = Array{Float64}(2,length(Ns),M,4,T)
  #for (methodidx,method) in enumerate([false,true])
  for (methodidx,method) in enumerate([pf,SQMC])
    println("Method: "*string(method))
    xhat, xreal = 0, 0
    @time for (i, N) in enumerate(Ns)
      println(string(N)*", ")
      rmslocal = Array{Float64,1}(M)
      for j = 1:M
        print(string(j)*",")
        #xreal, xhat = runTest(N,pf,debug,T,method)
        #xreal, xhat = runTest(N,pf,twoBeams2DSettings(method),T,debug)
        if method == pf
          xreal, xhat = runTest(N,pf,twoBeams2DSettings(),T,debug)
        else
          xreal, xhat = runTest(N,SQMC,twoBeams2DSettingsSQMC(),T,debug)
        end
        allxhats[methodidx,i,j,1:4,:] = xhat[[1,2,5,6],:]
        RMS[i,j,methodidx] = rms(xreal[1,:]-xhat[1,:])
      end
    end
    rmsMean[:,methodidx] = mean(RMS[:,:,methodidx],2)
    rmsVariance[:,methodidx] = std(RMS[:,:,methodidx],2)./sqrt(M)
    for (i, N) in enumerate(Ns)
      largeRMS[i,methodidx] = length(find(RMS[i,:,methodidx].>rmsMean[i,methodidx]+2*rmsVariance[i,methodidx]))
    end
  end
  rmsMean, rmsVariance, Ns, largeRMS, RMS, allxhats
end

function testPlot(rmsMean, rmsVariance, Ns)
  p = plot(x=Ns,y=rmsMean[:,1],c=:red)
  plot!(x=Ns,y=rmsMean[:,2],c=:blue)
  plot!(x=Ns,y=rmsMean[:,1]+rmsVariance[:,1]*2,c=:red)
  plot!(x=Ns,y=rmsMean[:,1]-rmsVariance[:,1]*2,c=:red)
  plot!(x=Ns,y=rmsMean[:,2]+rmsVariance[:,2]*2,c=:blue)
  plot!(x=Ns,y=rmsMean[:,2]-rmsVariance[:,2]*2,c=:blue)
  gui(p)
end
#
# function heatBoxPlot(plt, x, t, minmax, nbinsy=30)
#   if maximum(x)-minimum(x) > 0
#     heatmap!(plt, x=(t-1.5):1/length(x):(t-0.5),y=x',nbins=(1,nbinsy), ylims=tuple(minmax...))
#   end
#   return
# end
#
# function distPlotY(plt, x, t, minmax, color=:blue, nbinsy=30)
#   if maximum(x)-minimum(x) > 0
#     y, nbr = hist(x', nbinsy)
#     xval = nbr[:]#cumsum(nbr[:])
#     plot!(plt, [0, xval[:]]/xval[end], y, linecolor=color, ylims=tuple(minmax...))
#   end
#   return
# end
#
# """ `plotPoints(x, w, y, yhat, N, a, t, xreal, xhat, xOld, pdata)`
#
# To be called inside a particle filter, plots either particle density (`density=true`) or individual particles (`density=false`) \n
# To only plot the particle trajectories, set (`leftOnly=false`)\n
# Will plot all the real states in `pltxIdx` as well as the expected vs real measurements of `pltyIdx`.
#
# Arguments: \n
#   * `x`: `Array(M,N)`. The states for each patyicle where `M` number of states, `N` number of Particles
#   * `w`: `Array(N)`. logprobability of each particle
#   * `y`: `Array(R,T)`. All true outputs. `R` is number of outputs, `T` is total number of time steps (will only use index `t`)
#   * `yhat`: `Array(R,N)` The expected output per particle. `R` is number of outputs, `N` number of Particles
#   * `N`, Number of particles
#   * `a`, `Array(N)`, reorderng of particles (e.g. `1:N`)
#   * `t`, Current time step
#   * `xreal`: `Array(M,T)`. All true states. `R` is number of states, `T` is total number of time steps (will only use index `t`)
#   * `xhat`: Not used
#   * `xOld`: Same as `x`, but for previous time step, only used when `!density` to show states origins
#   * `pdata`: Persistant data for plotting. Set to void in first call and pdataOut on remaining \n
# Returns: `pdataOut`
# """
# function plotPoints(x, w, y, yhat, N, a, t, xreal, xhat, xOld, pdata)
#   immerse()
#   leftOnly = false
#   density = true
#   cols = leftOnly?1:2
#   grd = (r,c) -> (r-1)*cols+c
#   println("Surviving: "*string((N-length(setdiff(Set(1:N),Set(a))))/N))
#   vals = [x;yhat]
#   realVals = [xreal;y]
#   valsOld = [xOld;yhat]
#   pltxIdx = [1,2,5, 6]# 7 9]
#   pltyIdx = [1,2]
#   pltIdx = [pltxIdx; size(x,1)+pltyIdx]
#   println(pltIdx)
#   if pdata == Void
#     pdata = (subplot(layout=cols*ones(Int,length(pltIdx))), Array{Float64,2}(length(pltIdx),2))
#     gui(pdata[1])
#   end
#   p, minmax = pdata
#   println(size(vals[pltIdx,:],2))
#   minmax = [min(minmax[:,1], minimum(vals[pltIdx,:],2)) max(minmax[:,2], maximum(vals[pltIdx,:],2))]
#   #c = exp(w[:]-minimum(w))*3
#   c = exp(w[:])*5*N
#   for (i, val) in enumerate(pltIdx)
#     if !leftOnly
#       oldFig = p.plts[grd(i,2)].o[1]
#       newPlot = plot()
#     end
#     if density
#       #Plot the heatmap on the left plot
#       heatBoxPlot(p.plts[grd(i,1)], vals[val,:], t, minmax[i,:])
#     end
#     for j = 1:N
#       if !density
#         #Plot the line on the left plot
#         plot!(p.plts[grd(i,1)], [t-1.5,t-1], [valsOld[val,j], valsOld[val,j]], legend=false)
#         plot!(p.plts[grd(i,1)], [t-1,t-0.5], [valsOld[val,a[j]], vals[val,j]], legend=false)
#       end
#       if !leftOnly && !density
#         #Plot each of the dots on the right side
#         plot!(newPlot, [j,j], [vals[val,j]',vals[val,j]'], marker=(:circle, :red, c[j])  , legend=false)
#       end
#     end
#     if !leftOnly
#       if density
#         distPlotY(newPlot, vals[val,:] , t, minmax[i,:], :blue)
#         #Weight can not be multiplied with value! DUH!
#         #distPlotY(newPlot, vals[val,:].*w' , t, :red)
#       end
#       #Fix the bug with updating plot
#       p.plts[grd(i,2)] = copy(newPlot)
#       p.plts[grd(i,2)].o = (oldFig, p.plts[grd(i,2)].o[2])
#     end
#     #Plot Real State Here
#     plot!(p.plts[grd(i,1)], (t-1):t, [realVals[val,t], realVals[val,t]], legend=false, color=:black, linewidth=5)
#   end
#   gui(p)
#   print("Waiting for command. q to Quit, ^D to run all:\n")
#   if readline(STDIN) == "q\n"
#     error("Exit command")
#   end
#   (p, minmax)
# end