text:
"""
    mpccmodel_build_fn_from_defn(dimspec, defnnum, x, pr, ps)

Assemble `f`, `ce`, `ci`, and `F` Julia native functions from Symbolic
definition and returns both non-indexed and indexed variants. Generally no need
to call directly as this is called from `mpccmodel_load_defn_from_file()`.

**NOTE:** To understand how this works, one should also play with a few examples
of`build_function()` in Symbolics.jl.  The output from `build_function()` is a
2-element vector whereby the first element is the standard function adn the
second is a mutating version of the same function. Hence the formulation in the
code here.

# Arguments
- `dimspec::MPCCDimSpec`: The dimension specifications for the proram.
- `defnnum::MPCCDefinition`: The Symbolics.jl Num definitions of the
  expressions.
- `x`, `pr`, `ps`: Symbolics.jl Num variables for x, pr, ps.

"""
function mpccmodel_build_fn_from_defn(dimspec::MPCCDimSpec, defnnum::MPCCDefinition, x, pr, ps)

    @unpack n, q, l, me, mi = dimspec

    # Before building any functions, we need to setup the Num for Fq (product down columns)
    Fq = Vector{Num}(undef, q)
    for lp_q=1:q
        Fq[lp_q] = prod(defnnum.F[:, lp_q])
    end

    # Objective fn
    f_fn = build_function( [ defnnum.f ], x, pr, ps; expression=Val{false} )

    # Whole vector / array functions
    ce_fn = build_function( defnnum.ce, x, pr, ps; expression=Val{false} )
    ci_fn = build_function( defnnum.ci, x, pr, ps; expression=Val{false} )
    if ( l == 0 || q == 0 )
        # Workaround for wee bug in Symbolics when trying to compile empty matrix
        # it produces keech. So we hand it an empty vector instead.
        F_fn = build_function( Vector{Num}([]), x, pr, ps; expression=Val{false} )
    else
        F_fn = build_function( defnnum.F, x, pr, ps; expression=Val{false} )
    end
    Fq_fn = build_function( Fq, x, pr, ps; expression=Val{false} )

    # Indexed functions (need to separate std from mutating versions with [1] and [2])
    ce_i_fn = Vector{Tuple{Function, Function}}(undef, me)
    for lp_ce=1:me
        ce_i_fn[lp_ce] = build_function( [ defnnum.ce[lp_ce] ], x, pr, ps; expression=Val{false} )
    end
    ce_i_fn_std = [ ce_i_fn[lp_ce][1] for lp_ce=1:me ]
    ce_i_fn_mut = [ ce_i_fn[lp_ce][2] for lp_ce=1:me ]

    ci_i_fn = Vector{Tuple{Function, Function}}(undef, mi)
    for lp_ci=1:mi
        ci_i_fn[lp_ci] = build_function( [ defnnum.ci[lp_ci] ], x, pr, ps; expression=Val{false} )
    end
    ci_i_fn_std = [ ci_i_fn[lp_ci][1] for lp_ci=1:mi ]
    ci_i_fn_mut = [ ci_i_fn[lp_ci][2] for lp_ci=1:mi ]

    F_i_fn = Matrix{Tuple{Function, Function}}(undef, l, q)
    for lp_q=1:q
        for lp_l=1:l
            F_i_fn[lp_l, lp_q] = build_function( [ defnnum.F[lp_l, lp_q] ], x, pr, ps; expression=Val{false} )
        end
    end
    F_i_fn_std = [ F_i_fn[lp_l, lp_q][1] for lp_l=1:l, lp_q=1:q ]
    F_i_fn_mut = [ F_i_fn[lp_l, lp_q][2] for lp_l=1:l, lp_q=1:q ]

    Fq_i_fn = Vector{Tuple{Function, Function}}(undef, q)
    for lp_q=1:q
        Fq_i_fn[lp_q] = build_function( [ Fq[lp_q] ], x, pr, ps; expression=Val{false} )
    end
    Fq_i_fn_std = [ Fq_i_fn[lp_q][1] for lp_q=1:q ]
    Fq_i_fn_mut = [ Fq_i_fn[lp_q][2] for lp_q=1:q ]

    return MPCCFunctions(
            f_fn[1], f_fn[2],
            ce_fn[1], ce_fn[2],
            ce_i_fn_std, ce_i_fn_mut,
            ci_fn[1], ci_fn[2],
            ci_i_fn_std, ci_i_fn_mut,
            F_fn[1], F_fn[2],
            F_i_fn_std, F_i_fn_mut,
            Fq_fn[1], Fq_fn[2],
            Fq_i_fn_std, Fq_i_fn_mut
        )
end

"""
    mpccmodel_build_parameterisation(defns, t)

Build parameterisation from `defns` and Symbolics Num `t`, returns struct with
Julia native functions including derivatives.

"""
function mpccmodel_build_parameterisation(defns::Vector{MPCCParameterisationDefn}, t::Num)

    len_defns = length(defns)

    fns = Vector{MPCCParameterisationFunctions}(undef, len_defns)
    bf_results = Vector{Any}(undef, len_defns)
    bf_dt_results = Vector{Any}(undef, len_defns)
    sym_prdt = Vector{Vector{Num}}(undef, len_defns)
    for lp_i=1:len_defns
        # Calc derivatives
        sym_prdt[lp_i] = Symbolics.derivative(defns[lp_i].pr, t)

        # Compile them
        bf_results[lp_i] = Symbolics.build_function(defns[lp_i].pr, t; expression=Val{false})
        bf_dt_results[lp_i] = Symbolics.build_function(sym_prdt[lp_i], t; expression=Val{false})
        fns[lp_i] = MPCCParameterisationFunctions(
                bf_results[lp_i][1],
                bf_results[lp_i][2],
                bf_dt_results[lp_i][1],
                bf_dt_results[lp_i][2]
            )
    end

    return MPCCParameterisations(t, defns, fns)
end

function mpccmodel_load_defn_from_file(model_id::String)

    # TODO some sanity checking on model_id

    # Filename of the model we want to import, include() it
    src_filename = joinpath(@__DIR__, "models", "model_spec_$(model_id).jl")
    println("Loading model_spec_$(model_id).jl")
    include(src_filename)

    # We need to call two specific functions in the included file.
    # Need to setup these function variables from strings then we'll later use invokelatest
    src_fn_dimspec = "mm_spec_" * model_id * "_dimspec"
    src_fn_defn = "mm_spec_" * model_id * "_defn"
    src_fn_nzmask = "mm_spec_" * model_id * "_nzmask"
    src_fn_testvectors = "mm_spec_" * model_id * "_testvectors"
    src_fn_knownsols = "mm_spec_" * model_id * "_knownsols"
    src_fn_param_defns = "mm_spec_" * model_id * "_parameterisations"
    fn_dimspec = getfield(@__MODULE__, Symbol(src_fn_dimspec))
    fn_defn = getfield(@__MODULE__, Symbol(src_fn_defn))
    fn_nzmask = getfield(@__MODULE__, Symbol(src_fn_nzmask))
    fn_testvectors = getfield(@__MODULE__, Symbol(src_fn_testvectors))
    fn_knownsols = getfield(@__MODULE__, Symbol(src_fn_knownsols))
    fn_param_defns = getfield(@__MODULE__, Symbol(src_fn_param_defns))

    # We're getting the dimspec and setting up Symbolics.jl variables
    dimspec = Base.invokelatest(fn_dimspec)
    @variables x[1:dimspec.n] pr[1:dimspec.r] ps[1:dimspec.s] t

    # Grab the model definition, expressed in Symbolics.jl Num
    defn = Base.invokelatest(fn_defn, x, pr, ps)

    # Grab the param definitions
    param_defns = Base.invokelatest(fn_param_defns, t)

    # Build callable functions for the basic model: f, ce, ci, F
    fns = mpccmodel_build_fn_from_defn(dimspec, defn, x, pr, ps)

    # Build parameterisations
    parameterisations = mpccmodel_build_parameterisation(param_defns, t)

    nzmask = Base.invokelatest(fn_nzmask)
    testvectors = Base.invokelatest(fn_testvectors)
    knownsols = Base.invokelatest(fn_knownsols)

    return MPCCModelConfig(
        Symbolics.scalarize(x), Symbolics.scalarize(pr), Symbolics.scalarize(ps),
        dimspec,
        defn,
        fns,
        nzmask,
        testvectors,
        knownsols,
        parameterisations
    )
end

# Because JuMP is a bit shit, we need to construct scalar functions for it
# We fix the parameter set here too.
function mpccmodel_build_fixed_jump_fns(model::MPCCModelConfig, pr0::Vector{R}, ps0::Vector{Int64}) where {R <: Real}

    # Create suitable objective function for use in JuMP, fixing parameters
    jump_f(y::T...) where {T <: Real} = model.fns.f(collect(y), pr0, ps0)[1]  # This would otherwise return a vector!

    # Create suitable penalty objective function for use in JuMP, fixing parameters
    jump_f_pen = function(ρ::R, y::T...) where {R <: Real, T <: Real}
        x = collect(y)

        # Initial f value
        f_val = model.fns.f(collect(y), pr0, ps0)[1]

        # Add penalty
        dot_prod = zero(T)
        for lp_q=1:model.dimspec.q
            cnstr_prod = one(T)
            for lp_l=1:model.dimspec.l
                cnstr_prod *= abs(model.fns.F_i[lp_l, lp_q](collect(y), pr0, ps0)[1])
            end
            dot_prod += cnstr_prod
        end

        return (f_val + ρ*dot_prod)
    end

    # Create suitable constraint functions for use in JuMP, fixing parameters
    jump_ce = Vector{Function}(undef, model.dimspec.me)
    for lp_me=1:model.dimspec.me
        jump_ce[lp_me] = function (y::T...) where {T <: Real}
            return model.fns.ce_i[lp_me](collect(y), pr0, ps0)[1]
        end
    end

    # Create suitable constraint functions for use in JuMP, fixing parameters
    jump_ci = Vector{Function}(undef, model.dimspec.mi)
    for lp_mi=1:model.dimspec.mi
        jump_ci[lp_mi] = function (y::T...) where {T <: Real}
            return model.fns.ci_i[lp_mi](collect(y), pr0, ps0)[1]
        end
    end

    # Create suitable constraint functions for use in JuMP, fixing parameters
    jump_F = Matrix{Function}(undef, model.dimspec.l, model.dimspec.q)
    for lp_q=1:model.dimspec.q
        for lp_l=1:model.dimspec.l
            jump_F[lp_l, lp_q] = function (y::T...) where {T <: Real}
                return model.fns.F_i[lp_l, lp_q](collect(y), pr0, ps0)[1]
            end
        end
    end

    return MPCCJumpFixedFunctions(jump_f, jump_f_pen, jump_ce, jump_ci, jump_F)
end


text:
# Types

struct AWAP <: SDMDataSource end

struct Temperature{X} end
struct VapourPressure{X} end
struct Solar{X} end
struct Rainfall{X} end

struct H09 end
struct H15 end
struct MinAve end
struct MaxAve end

const AWAP_DATEFORMAT = DateFormat("yyyymmdd")

# Interface methods

download_raster(T::Type{AWAP}, layers::Tuple=rasterlayers(T); kwargs...) =
    map(l -> download_raster(T, l; kwargs...), layers)
function download_raster(T::Type{AWAP}, layer::Type; dates)
    dates = _date_sequence(dates, Day(1))
    mkpath(rasterpath(T, layer))
    raster_paths = String[]
    for d in dates
        raster_path = rasterpath(T, layer, d)
        if !isfile(raster_path)
            zip_path = zippath(T, layer, d)
            _maybe_download(zipurl(T, layer, d), zip_path)
            run(`uncompress $zip_path -f`)
        end
        push!(raster_paths, raster_path)
    end
    return raster_paths
end

rasterlayers(::Type{<:AWAP}) =
    (Solar, Rainfall, VapourPressure{H09}, VapourPressure{H15},
     Temperature{MinAve}, Temperature{MaxAve})

rasterpath(T::Type{AWAP}) = joinpath(rasterpath(), "AWAP")
rasterpath(T::Type{AWAP}, layer) = joinpath(rasterpath(T), _pathsegments(layer)[1:2]...)
rasterpath(T::Type{AWAP}, layer, date::Dates.AbstractTime) =
    joinpath(rasterpath(T, layer), rastername(T, layer, date))

rastername(T::Type{AWAP}, layer, date::Dates.AbstractTime) =
    joinpath(_date2string(T, date) * ".grid")

function zipurl(T::Type{AWAP}, layer, date)
    s = _pathsegments(layer)
    d = _date2string(T, date)
    # The actual zip name has the date twice, which is weird.
    # So we download in to a different name as there no output
    # name flages for `uncompress`. It's ancient.
    joinpath("http://www.bom.gov.au/web03/ncc/www/awap", s..., "grid/0.05/history/nat/$d$d.grid.Z")
end

zipname(T::Type{AWAP}, layer, date) = _date2string(T, date) * ".grid.Z"

zippath(T::Type{AWAP}, layer, date) = joinpath(rasterpath(T, layer), zipname(T, layer, date))

# Utilitiy methods

_dateformat(::Type{AWAP}) = DateFormat("yyyymmdd")

_pathsegments(::Type{<:Solar}) = "solar", "solarave", "daily"
_pathsegments(::Type{<:Rainfall}) = "rainfall", "totals", "daily"
_pathsegments(::Type{VapourPressure{H09}}) = "vprp", "vprph09", "daily"
_pathsegments(::Type{VapourPressure{H15}}) = "vprp", "vprph15", "daily"
_pathsegments(::Type{Temperature{MinAve}}) = "temperature", "minave", "daily"
_pathsegments(::Type{Temperature{MaxAve}}) = "temperature", "maxave", "daily"

#=
Add ndvi monthly?
ndvi,ndviave,month
=#

ext: jl


text:
#############################################################################
#Copyright (c) 2017 Charles Le Losq
#
#The MIT License (MIT)
#
#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.
#
#############################################################################

VERSION >= v"0.5.0" && __precompile__()

module PeriodicTable

using DataFrame

include("mass.jl")

#From mass.jl
export mass

end # module


text:
using MixedModels: fixef!, stderror!
using MixedModels: getθ!, setθ!, updateL!
using LinearAlgebra
using Statistics

"""
    permutation([rng::AbstractRNG,] nsamp::Integer, m::LinearMixedModel;
                use_threads::Bool=false,
                β=zeros(length(coef(morig))),
                residual_method=:signflip,
                blup_method=ranef)

Perform `nsamp` nonparametric bootstrap replication fits of `m`, returning a `MixedModelBootstrap`.

The default random number generator is `Random.GLOBAL_RNG`.

`GeneralizedLinearMixedModel` is currently unsupported.

# Named Arguments
`use_threads` determines whether or not to use thread-based parallelism.

!!! note
    Note that `use_threads=true` may not offer a performance boost and may even
    decrease peformance if multithreaded linear algebra (BLAS) routines are available.
    In this case, threads at the level of the linear algebra may already occupy all
    processors/processor cores. There are plans to provide better support in coordinating
    Julia- and BLAS-level threads in the future.

!!! warning
    The PRNG shared between threads is locked using [`Threads.SpinLock`](@ref), which
    should not be used recursively. Do not wrap `permutation` in an outer `SpinLock`.

`hide_progress` can be used to disable the progress bar. Note that the progress
bar is automatically disabled for non-interactive (i.e. logging) contexts.

Permutation at the level of residuals can be accomplished either via sign
flipping (`residual_method=:signflip`) or via classical
permutation/shuffling (`residual_method=:shuffle`).

`blup_method` provides options for how/which group-level effects are passed for permutation.
The default `ranef` uses the shrunken conditional modes / BLUPs. Unshrunken estimates from
ordinary least squares (OLS) can be used with `olsranef`. There is no shrinkage of the
group-level estimates with this approach, which means singular estimates can be avoided.
However, if the design matrix for the random effects is rank deficient (e.g., through the use
of `MixedModels.fulldummy` or missing cells in the data), then this method will fail.
See [`olsranef`](@ref) and `MixedModels.ranef` for more information.

Generally, permutations are used to test a particular (null) hypothesis. This
hypothesis is specified via by setting `β` argument to match the hypothesis. For
example, the null hypothesis that the first coefficient is zero would expressed as

```julialang
julia> hypothesis = coef(model);
julia> hypothesis[1] = 0.0;
```
!!! note
    The permutation (test) generates samples from H0, from which
    it is possible to compute p-values. The bootstrap typically generates
    samples from H1, which are convenient for generating coverage/confidence
    intervals. Of course, confidence intervals and p-values are duals of each
    other, so it is possible to convert from one to the other.

# Method

The method implemented here is based on the approach given in:

ter Braak C.J.F. (1992) Permutation Versus Bootstrap Significance Tests in
Multiple Regression and Anova. In: Jöckel KH., Rothe G., Sendler W. (eds)
Bootstrapping and Related Techniques. Lecture Notes in Economics and Mathematical
Systems, vol 376. Springer, Berlin, Heidelberg.
https://doi.org/10.1007/978-3-642-48850-4_10

and

Winkler, A. M., Ridgway, G. R., Webster, M. A., Smith, S. M., & Nichols, T. E. (2014).
Permutation inference for the general linear model. NeuroImage, 92, 381–397.
https://doi.org/10.1016/j.neuroimage.2014.01.060

!!! warning
    This method has serious limitations for singular models because sign-flipping a zero
    is not an effective randomization technique.
"""
function permutation(
    rng::AbstractRNG,
    n::Integer,
    morig::LinearMixedModel{T};
    use_threads::Bool=false,
    hide_progress=false,
    β::AbstractVector{T}=zeros(T, length(coef(morig))),
    residual_method=:signflip,
    blup_method=ranef_scaled,
) where {T}
    # XXX instead of straight zeros,
    #     should we use 1-0s for intercept only?
    βsc, θsc = similar(morig.β), similar(morig.θ)
    p, k = length(βsc), length(θsc)
    m = deepcopy(morig)

    β_names = (Symbol.(fixefnames(morig))..., )
    rank = length(β_names)

    blups,scalings = blup_method(morig)

    resids = residuals(morig)#, blups)
    reterms = morig.reterms

    # we need arrays of these for in-place operations to work across threads
    m_threads = [m]
    βsc_threads = [βsc]
    θsc_threads = [θsc]

    if use_threads
        Threads.resize_nthreads!(m_threads)
        Threads.resize_nthreads!(βsc_threads)
        Threads.resize_nthreads!(θsc_threads)
    end
    # we use locks to guarantee thread-safety, but there might be better ways to do this for some RNGs
    # see https://docs.julialang.org/en/v1.3/manual/parallel-computing/#Side-effects-and-mutable-function-arguments-1
    # see https://docs.julialang.org/en/v1/stdlib/Future/index.html
    rnglock = Threads.SpinLock()
    samp = replicate(n; use_threads=use_threads, hide_progress=hide_progress) do
        tidx = use_threads ? Threads.threadid() : 1
        model = m_threads[tidx]
        local βsc = βsc_threads[tidx]
        local θsc = θsc_threads[tidx]
        lock(rnglock)
        model = permute!(rng, model; β=β, blups=blups, resids=resids,
                         residual_method=residual_method, scalings=scalings)
        unlock(rnglock)
        refit!(model)
        (
         objective = model.objective,
         σ = model.σ,
         β = NamedTuple{β_names}(fixef!(βsc, model)),
         se = SVector{p,T}(stderror!(βsc, model)),
         θ = SVector{k,T}(getθ!(θsc, model)),
        )
    end
    MixedModelPermutation(
        samp,
        deepcopy(morig.λ),
        getfield.(morig.reterms, :inds),
        copy(morig.optsum.lowerbd),
        NamedTuple{Symbol.(fnames(morig))}(map(t -> (t.cnames...,), morig.reterms)),
    )
end

function ranef_scaled(morig)
    blups = ranef(morig)
    resids = residuals(morig)
    scalings = inflation_factor(morig, blups, resids)
    return blups,scalings
end

function permutation(nsamp::Integer, m::LinearMixedModel, args...; kwargs...)
    return permutation(Random.GLOBAL_RNG, nsamp, m, args...; kwargs...)
end

function permutation(rng::AbstractRNG, n::Integer,
                                morig::GeneralizedLinearMixedModel;
                                use_threads::Bool=false) where {T}

    throw(ArgumentError("GLMM support is not yet implemented"))
end

permute!(model::LinearMixedModel, args...; kwargs...) =
    permute!(Random.GLOBAL_RNG, model, args...; kwargs...)

"""
    permute!([rng::AbstractRNG,] model::LinearMixedModel;
             β=zeros(length(coef(model))),
             blups=ranef(model),
             resids=residuals(model,blups),
             residual_method=:signflip,
             scalings=inflation_factor(model))

Simulate and install a new response via permutation of the residuals
at the observational level and sign-flipping of the conditional modes at group level.

Generally, permutations are used to test a particular (null) hypothesis. This
hypothesis is specified via by setting `β` argument to match the hypothesis. For
example, the null hypothesis that the first coefficient is zero would expressed as

```julialang
julia> hypothesis = coef(model);
julia> hypothesis[1] = 0.0;
```

Permutation at the level of residuals can be accomplished either via sign
flipping (`residual_method=:signflip`) or via classical
permutation/shuffling (`residual_method=:shuffle`).

Sign-flipped permutation of the residuals is similar to permuting the
(fixed-effects) design matrix; shuffling the residuals is the same as permuting the
(fixed-effects) design matrix. Sign-flipping the random effects
preserves the correlation structure of the random effects, while also being equivalent to
permutation via swapped labels for categorical variables.

!!! warning
    This method has serious limitations for singular models because sign-flipping a zero
    is not an effective randomization technique.

Optionally, instead of using the shrunken random effects from `ranef`, within-group OLS
estimates can be computed and used instead with [`olsranef`](@ref). There is no shrinkage
of the group-level estimates with this approach, which means singular estimates can be
avoided. However, if the design matrix for the random effects is rank deficient (e.g.,
through the use of `MixedModels.fulldummy` or missing cells in the data), then this method
will fail.

In addition to the permutation step, there is also an inflation step. Due to the
shrinkage associated with the random effects in a mixed model, the variance of the
conditional modes / BLUPs / random intercepts and slopes is less than the variance
estimated by the model and displayed in the model summary or via `MixedModels.VarCorr`.
This shrinkage also impacts the observational level residuals. To compensate for this,
the resampled residuals and groups are scale-inflated so that their standard deviation
matches that of the estimates in the original model. The default inflation factor is
computed using [`inflation_factor`](@ref) on the model.

See also [`permutation`](@ref), [`nonparametricbootstrap`](@ref) and [`resample!`](@ref).

The functions `coef` and `coefnames` from `MixedModels` may also be useful.

# Reference
The method implemented here is based on the approach given in:

ter Braak C.J.F. (1992) Permutation Versus Bootstrap Significance Tests in
Multiple Regression and Anova. In: Jöckel KH., Rothe G., Sendler W. (eds)
Bootstrapping and Related Techniques. Lecture Notes in Economics and Mathematical
Systems, vol 376. Springer, Berlin, Heidelberg.
https://doi.org/10.1007/978-3-642-48850-4_10

and

Winkler, A. M., Ridgway, G. R., Webster, M. A., Smith, S. M., & Nichols, T. E. (2014).
Permutation inference for the general linear model. NeuroImage, 92, 381–397.
https://doi.org/10.1016/j.neuroimage.2014.01.060
"""
function permute!(rng::AbstractRNG, model::LinearMixedModel{T};
                  β::AbstractVector{T}=zeros(T, length(coef(model))),
                  blups=ranef(model),
                  resids=residuals(model,blups),
                  residual_method=:signflip,
                  scalings=inflation_factor(model)) where {T}

    reterms = model.reterms
    y = response(model) # we are now modifying the model
    copy!(y, resids)

    # inflate these to be on the same scale as the empirical variation instead of the MLE
    y .*= last(scalings)

    if residual_method == :shuffle
        shuffle!(rng, y)
    elseif  residual_method == :signflip
        y .*= rand(rng, (-1,1), length(y))
    else
        throw(ArgumentError("Invalid: residual permutation method: $(residual_method)"))
    end

    for (inflation, re, trm) in zip(scalings, blups, reterms)
        npreds, ngrps = size(re)
        # sign flipping
        newre = re * diagm(rand(rng, (-1,1), ngrps))

        # this just multiplies the Z matrices by the BLUPs
        # and add that to y
        mul!(y, trm, lmul!(inflation, newre), one(T), one(T))
        # XXX inflation is resampling invariant -- should we move it out?
    end

    mul!(y, model.X, β, one(T), one(T))

    # mark model as unfitted
    model.optsum.feval = -1

    return model
end

permutationtest(perm::MixedModelPermutation, model::LinearMixedModel) = permutationtest(perm::MixedModelPermutation, model, :twosided)

"""
    permutationtest(perm::MixedModelPermutation, model, type=:greater)

Perform a permutation using the already computed permutation and given the observed values.

The `type` parameter specifies use of a two-sided test (`:twosided`) or the directionality of a one-sided test
(either `:lesser` or `:greater`, depending on the hypothesized difference to the null hypothesis).

See also [`permutation`](@ref).

To account for finite permutations, we implemented the conservative method from Phipson & Smyth 2010:
 Permutation P-values Should Never Be Zero:Calculating Exact P-values When Permutations Are Randomly Drawn
 http://www.statsci.org/webguide/smyth/pubs/permp.pdf

"""
function permutationtest(perm::MixedModelPermutation, model; type::Symbol=:twosided,β::AbstractVector=zeros(length(coef(model))), statistic=:z)
    #@warn """This method is known not to be fully correct.
    #         The interface for this functionality will likely change drastically in the near future."""
    # removed due to distributed run

    if type == :greater || type  == :twosided
        comp = >=
    elseif type == :lesser
        comp = <=
    else
        throw(ArgumentError("Comparison type $(type) unsupported"))
    end
    if statistic == :z
        x = coeftable(model)
        ests = Dict(Symbol(k) => v for (k,v) in zip(coefnames(model), x.cols[x.teststatcol]))
    elseif statistic == :β
        ests = Dict(Symbol(k) => v for (k,v) in zip(coefnames(model), coef(model)))
    else
        error("statistic not implemented yet")
    end

    perms = columntable(perm.coefpvalues)

    dd = Dict{Symbol, Vector}()

    for (ix,k) in enumerate(Symbol.(coefnames(model)))
        dd[k] = perms[statistic][perms.coefname .== k]

        push!(dd[k],ests[k]) # simplest approximation to ensure p is never 0 (impossible for permutation test)
        if type == :twosided
            # in case of testing the betas, H0 might be not β==0, therefore we have to remove it here first before we can abs
            # the "z's" are already symmetric around 0 regardless of hypothesis.
            if statistic == :β
                #println(β[ix])
                dd[k]  .= dd[k]  .- β[ix]
                ests[k] = ests[k] - β[ix]
            end

              dd[k]  .= abs.(dd[k])
              ests[k] = abs(ests[k])
        end

    end

    # short way to calculate:
    # b = sum.(abs.(permDist).>=abs.(testValue)); (twosided)
    # Includes the conservative correction for approximate permutation tests
    # p_t = (b+1)/(nperm+1);

    # (with comp being <=) Note that sum(<=(ests),v) does the same as  sum(v .<=ests) (thus "reversed" arguments in the first bracket)
    results = (; (k => (1+sum(comp(ests[k]),v))/(1+length(v)) for (k,v) in dd)...)

    return results
end


text:
import brainflow

# specify logging library to use
brainflow.enable_dev_brainflow_logger(Integer(brainflow.BOARD_CONTROLLER))

params = brainflow.BrainFlowInputParams()
board_shim = brainflow.BoardShim(Integer(brainflow.SYNTHETIC_BOARD), params)
sampling_rate = brainflow.get_sampling_rate(Integer(brainflow.SYNTHETIC_BOARD))
nfft = brainflow.get_nearest_power_of_two(sampling_rate)

brainflow.prepare_session(board_shim)
brainflow.start_stream(board_shim)
sleep(5)
brainflow.stop_stream(board_shim)
data = brainflow.get_board_data(board_shim)
brainflow.release_session(board_shim)

eeg_channels = brainflow.get_eeg_channels(Integer(brainflow.SYNTHETIC_BOARD))

bands = brainflow.get_avg_band_powers(data, eeg_channels, sampling_rate, true)
println(bands[1])
println(bands[2])

max_stars_repo_name: Cognixion-inc/brainflow

max_issues_repo_name: Cognixion-inc/brainflow

max_forks_repo_name: Cognixion-inc/brainflow


text:
function main()
    gross = 0
    a = 1
    b = 2

    while a < 4e6
        if a % 2 == 0
            gross += a
        end
        c = a + b
        a = b
        b = c
    end

    println(gross)
end

main()

max_stars_repo_name: shanedrabing/project-euler

max_issues_repo_name: shanedrabing/project-euler

max_forks_repo_name: shanedrabing/project-euler


text:
# function make_test_entry_results(switches :: AbstractEntrySwitches{F1}) where F1
# 	J = n_types(switches);
# 	nc = n_colleges(switches);
# 	eVal_jV = collect(range(-0.5, 1.5, length = J));
#     probEnter_jcM = range(0.1, 0.2, length = J) * range(1.0, 1.2, length = nc)';
# 	enrollV = college_enrollment(switches, probEnter_jcM);
# 	return EntryResults(switches, probEnter_jcM, eVal_jV, enrollV);
# end

## -------------  Construction

"""
	$(SIGNATURES)

Object that holds results of sequential entry protocol. If there is only one location, all students are considered "local."
"""
function EntryResults(switches :: EntryDecisionSwitches{F1}) where F1
	J = n_types(switches);
	nc = n_colleges(switches);
	nl = n_locations(switches);
	return EntryResults(switches,
		zeros(F1, J, nl, nc), zeros(F1, J, nl, nc), zeros(F1, J, nl, nc),
		zeros(F1, J, nl),
		zeros(F1, nc, nl), zeros(F1, nc, nl), zeros(F1, nc, nl))
end

"""
	$(SIGNATURES)

Validate `EntryResults`.
For subsetted `EntryResults`, do not `validateFracLocal`.
"""
function validate_er(er :: EntryResults{F1}; validateFracLocal :: Bool = true) where F1
	isValid = true
	isValid = isValid  &&  check_prob_array(er.fracEnter_jlcM);
	isValid = isValid  &&  check_prob_array(er.fracLocal_jlcM);
	isValid = isValid  &&  all(er.fracLocal_jlcM .< er.fracEnter_jlcM .+ 1e-5);
	isValid = isValid  &&  check_prob_array(sum(er.fracEnter_jlcM, dims = 3));
	if any(er.enrollLocal_clM .> er.enroll_clM .+ 1e-5)
		isValid = false;
		@warn "$er: Local enrollment > total enrollment"
	end

	nc = size(er.enroll_clM, 1);
	if any(er.enroll_clM .< er.enrollBest_clM .- 0.001)
		@warn "Total enrollment should be larger than best enrollment"
		isValid = false;
	end
	if !isapprox(er.enroll_clM[nc,:], er.enrollBest_clM[nc,:])
		@warn "For top college, total enrollment should equal best enrollment"
		isValid = false;
	end
	if any(map((x,y) -> x > y + 0.001,  er.fracEnterBest_jlcM, er.fracEnter_jlcM))
		maxGap, maxIdx = findmax(er.fracEnterBest_jlcM .- er.fracEnter_jlcM);
		j = maxIdx[1];
		@warn """
			More students in best college than total.
			Max gap: $maxGap for index $maxIdx
			fracEnterBest: $(er.fracEnterBest_jlcM[j,:,:])
			fracEnter:     $(er.fracEnter_jlcM[j,:,:])
			"""
		isValid = false;
	end
	if !isapprox(er.fracEnterBest_jlcM[:,:,nc], er.fracEnter_jlcM[:,:,nc]; atol = 1e-7)
		maxGap = maximum(abs.(er.fracEnterBest_jlcM[:,:,nc] .- er.fracEnter_jlcM[:,:,nc]));
		@warn """
			For top college, total entry should equal best entry.
			Max gap: $maxGap
			"""
		# isValid = false;  # enable again
	end

	# Computing fracLocal across colleges and across types should give the same answer
	if validateFracLocal
		fracLocal = frac_local(er);
		fracLocal2 = sum(frac_local_c(er) .* enrollment_c(er)) / sum(enrollment_c(er));
		entryMass_jV = entry_probs_j(er, :all) .* type_mass_j(er);
		fracLocal3 = sum(frac_local_j(er) .* entryMass_jV) / sum(entryMass_jV);
		if !isapprox(fracLocal, fracLocal2)  ||  !isapprox(fracLocal3, fracLocal)
			@warn "Fraction local not consistent: $fracLocal, $fracLocal2, $fracLocal3"
			isValid = false;
		end
	end
	return isValid
end

Base.show(io :: IO, er :: AbstractEntryResults{F1}) where F1 =
	print(io, typeof(er))

# Cannot make EntryResults without actually solving. At least its hard to maintain consistency.
function make_test_entry_results(switches :: EntryDecisionSwitches{F1};
	typicalValue :: F1 = 1.0) where F1
    rng = MersenneTwister(34);

	st = make_test_symbol_table();
	nc = n_colleges(switches);
	J = n_types(switches);
	nl = n_locations(switches);

    admissionS = make_test_admissions_cutoff(nc);
    objId = ObjectId(:entryOneStep);
    entryS = init_entry_decision(objId, switches, st);
	admProbFct = make_test_admprob_fct_logistic(nc);

	vWork_jV, vCollege_jcM = values_for_test(rng, J, nc, nl;
		typicalValue = typicalValue);
    hsGpaPctV = collect(range(0.1, 0.9, length = J));
    rank_jV = vcat(2 : 2 : J, 1 : 2 : J);

    er = entry_decisions(entryS, admissionS, admProbFct,
        vWork_jV, vCollege_jcM, hsGpaPctV, rank_jV);
	# er = EntryResults(switches);
	# for j = 1 : J
	# 	for l = 1 : nl
	# 		er.eVal_jlM[j, l] = 0.5 * j + 0.6 * l;
	# 		for ic = 1 : nc
	# 			er.fracEnter_jlcM[j, l, ic] = 0.01 * J + 0.02 * l + 0.015 * ic;
	# 			er.fracLocal_jlcM[j, l, ic] = 0.5 * er.fracEnter_jlcM[j, l, ic];
	# 		end
	# 	end
	# end
	# for l = 1 : nl
	# 	for ic = 1 : nc
	# 		er.enroll_clM[ic, l] = 0.3 * ic + 0.2 * l;
	# 		er.enrollLocal_clM[ic, l] = 0.3 * er.enroll_clM[ic, l];
	# 	end
	# end
	@assert validate_er(er)
	return er
end

## ----------  Helpers

# Take the mean of a variable by (j,l,c) across locations, weighted by typeMass(j,l).
function mean_over_locations(e :: EntryResults{F1}, x_jlcM :: Array{F1, 3}) where F1
	J, nl, nc = size(x_jlcM);
	x_jcM = zeros(F1, J, nc);
	for ic = 1 : nc
		x_jcM[:, ic] = sum(x_jlcM[:,:,ic] .* type_mass_jl(e),  dims = 2) ./
			type_mass_j(e);
	end
	return x_jcM
end

function mean_over_locations(e :: EntryResults{F1}, x_jlM :: Matrix{F1}) where F1
	x_jV = sum(x_jlM .* type_mass_jl(e),  dims = 2) ./
		type_mass_j(e);
	return x_jV
end

"""
	$(SIGNATURES)

Return a copy of the `EntryDecision` object for a subset of types.
College enrollments are not updated. There is not enough info to compute local enrollments.
"""
function subset_types(er :: EntryResults{F1}, idxV :: AbstractVector) where F1
	newSwitches = deepcopy(er.switches);
	subset_types!(newSwitches, idxV);
	erOut = EntryResults(newSwitches, er.fracEnter_jlcM[idxV,:,:],
		er.fracLocal_jlcM[idxV,:,:], er.fracEnterBest_jlcM[idxV,:,:],
		er.eVal_jlM[idxV,:],
		er.enroll_clM, er.enrollLocal_clM, er.enrollBest_clM);
	@assert validate_er(erOut; validateFracLocal = false)
	return erOut
end

## ----------  College enrollment

"""
	$(SIGNATURES)

College enrollments. By location (if any).
"""
function enrollment_cl(e :: AbstractEntryResults{F1}, univ :: Symbol = :all) where F1
	if univ == :all
		return e.enroll_clM;
	elseif univ == :local
		return e.enrollLocal_clM;
	elseif univ == :nonlocal
		return enrollment_cl(e, :all) .- enrollment_cl(e, :local);
	else
		error("Invalid $univ");
	end
end

enrollment_cl(e :: EntryResults{F1}, ic :: Integer,
	univ :: Symbol = :all) where F1 =
	enrollment_cl(e, univ)[ic,:];

enrollment_cl(e :: EntryResults{F1}, ic :: Integer, l,
	univ = :all) where F1 =
	enrollment_cl(e, univ)[ic,l];

"""
	$(SIGNATURES)

Total enrollment by college; across all locations.
"""
enrollment_c(e :: AbstractEntryResults{F1}, univ :: Symbol = :all) where F1 =
	vec(sum(enrollment_cl(e, univ), dims = 2));

# Mass of students in each college for who this is the best college
enroll_best_c(e :: AbstractEntryResults{F1}) where F1 =
	vec(sum(e.enrollBest_clM, dims = 2));

"""
	$(SIGNATURES)

Fraction of students in each college type for who this is the best available college (across all locations). For the top college, this is 1 by construction.
"""
frac_best_c(e :: AbstractEntryResults{F1}) where F1 =
	enroll_best_c(e) ./ enrollment_c(e);

## ----------  Fraction local

"""
	$(SIGNATURES)

Fraction of students who attend local colleges.
"""
function frac_local(e :: EntryResults{F1}) where F1
	if n_locations(e) == 1
		fracLocal = one(F1);
	else
		fracLocal = sum(frac_local_c(e) .* enrollment_c(e)) /
			sum(enrollment_c(e));
		# Prevent rounding errors
		fracLocal = min(fracLocal, one(F1));
	end
	return fracLocal
end

"""
	$(SIGNATURES)

Fraction of students of each type who attend local colleges.
"""
function frac_local_j(e :: EntryResults{F1}) where F1
	if n_locations(e) == 1
		fracLocalV = ones(F1, n_types(e));
	else
		fracLocalV = entry_probs_j(e, :local) ./ entry_probs_j(e, :all);
		fracLocalV = min.(fracLocalV, ones(F1));
	end
	return fracLocalV
end

"""
	$(SIGNATURES)

Fraction of students at each college who are local. Across all locations.
"""
function frac_local_c(e :: EntryResults{F1}) where F1
	if n_locations(e) == 1
		fracLocalV = ones(F1, n_colleges(e));
	else
		fracLocalV = enrollment_c(e, :local) ./ enrollment_c(e, :all);
		fracLocalV = min.(fracLocalV, ones(F1));
	end
	return fracLocalV
end

## --------------  Entry probabilities of types

"""
	$(SIGNATURES)

Entry probabilities by college for each type, location.
"""
function entry_probs_jlc(e :: EntryResults{F1}, univ :: Symbol = :all) where F1
	if univ == :local
		prob_jlcM = e.fracLocal_jlcM;
	elseif univ == :nonlocal
		prob_jlcM = e.fracEnter_jlcM .- e.fracLocal_jlcM;
	elseif univ == :all
		prob_jlcM = e.fracEnter_jlcM;
	else
		error("Invalid: $univ")
	end
end

"""
	$(SIGNATURES)

Fraction of each type entering each college; across all locations.
"""
function entry_probs_jc(e :: EntryResults{F1}, univ :: Symbol = :all) where F1
	mean_over_locations(e, entry_probs_jlc(e, univ));
end

"""
	$(SIGNATURES)

Entry probabilities by type, across all colleges.
"""
function entry_probs_j(er :: AbstractEntryResults{F1},
	univ :: Symbol = :all) where F1
	return vec(sum(entry_probs_jc(er, univ), dims = 2));
end

## ----------------  Expected values

"""
	$(SIGNATURES)

Expected value at college entry stage, by type.
"""
expected_values_jl(e :: EntryResults{F1}) where F1 = e.eVal_jlM;

"""
	$(SIGNATURES)

Expected values at college entry by type; averaged across locations.
"""
expected_values_j(e :: EntryResults{F1}) where F1 =
	mean_over_locations(e, expected_values_jl(e));

expected_values_j(e :: EntryResults{F1}, l :: Integer) where F1 =
	e.eVal_jlM[:, l];

"""
	$(SIGNATURES)

Returns `Vector{Bool}` that indicates which colleges are full.
"""
colleges_full(e :: AbstractEntryResults{F1}) where F1 =
    enrollment_cl(e) .>= capacities(e);

## -----------  Modify

"""
	$(SIGNATURES)

Scale entry probs to bound within `min_entry_prob` and `max_entry_prob`.
"""
function scale_entry_probs!(er :: AbstractEntryResults{F1}) where F1
    minEntryProb = min_entry_prob(er.switches);
	maxEntryProb = max_entry_prob(er.switches);
	# This creates the risk that `fracEnterBest` may no longer be consistent with `fracEnter`. To avoid this, `fracEnterBest` is scaled the complicated way.
	# The `min.(one, ...)` deals with numerical inaccuracies.
	bestToAll_jlcM = min.(one(F1),
		er.fracEnterBest_jlcM ./ max.(minEntryProb, er.fracEnter_jlcM));
	# Reach into object to ensure that we don't get a copy
	scale_entry_probs!(er.fracEnter_jlcM, minEntryProb, maxEntryProb);
	scale_entry_probs!(er.fracLocal_jlcM, minEntryProb, maxEntryProb);
	er.fracEnterBest_jlcM .= er.fracEnter_jlcM .* bestToAll_jlcM;
	# Imposing this directly avoids errors during validation
	nc = size(er.fracEnter_jlcM, 3);
	er.fracEnterBest_jlcM[:,:,nc] .= er.fracEnter_jlcM[:,:,nc];

	@assert all_less(er.fracEnterBest_jlcM, er.fracEnter_jlcM; atol = 0.0001)
	return nothing
end

"""
	$(SIGNATURES)

Scale entry probs to match given total entry probs by type.
Respect the fraction of each type that goes to each (l, c). Also adjust `fracLocal_jlcM` to ensure that the ratio of local to total entry is unchanged.

`EntryResults` are no longer internally consistent after scaling. Only meant for testing / fixing poor outcomes during an optimization.
"""
function fix_type_entry_probs!(er :: AbstractEntryResults{F1}, typeTotalV :: AbstractVector{F1}) where F1 <: Real

    @assert all_at_most(typeTotalV, one(F1))
	@assert length(typeTotalV) == n_types(er)

	# Ratio local/total entry
	fracEnter_jV = entry_probs_j(er, :all);

	J = n_types(er);
	for j = 1 : J
		er.fracEnter_jlcM[j, :, :] .*= (typeTotalV[j] / fracEnter_jV[j]);
		er.fracLocal_jlcM[j, :, :] .*= (typeTotalV[j] / fracEnter_jV[j]);
	end

	return nothing
end

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


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

using Markdown
using InteractiveUtils

# ╔═╡ 674905ed-75a1-457c-a866-279cdb76d243
#hide
begin
#for sharing environments between notebooks
import Pkg
Pkg.activate(Base.current_project())
end

# ╔═╡ e3cd4e1f-c7ac-45fd-af9c-1327c0f7a426
#hide
using PlutoUI

# ╔═╡ 7a0b47e0-0231-4d3f-b7b8-f6dda83f120d
#export
using OrderedCollections

# ╔═╡ 7831e355-bd3b-4333-8200-b8938277af3e
using PlutoTest

# ╔═╡ 2c396254-932f-4431-8f94-08dd2766aa8e
#hide
include("../src/Export.jl")

# ╔═╡ 17a13720-29ba-11ec-0ba4-83c9ee9ac582
#export
include("../src/Documenter.jl")

# ╔═╡ 3e8524c9-1b00-4e30-bb47-524991f3d0f7
#export
import YAML

# ╔═╡ 2ac37248-81d0-4c8b-99c7-4dd9f400a955
#hide
TableOfContents()

# ╔═╡ 1f7d1cb4-7fb9-4f40-9350-70ef40bca38f
md"## build"

# ╔═╡ 1c5d5f39-0349-454b-aa22-5feda1287b9d
#export
"""
> build(nbsdir)--> extracts markdown items and builds documentation from the notebooks in the given directory. Also, creates creates a readme file with content in the project root.
"""
function build(nbsdir)
    Documenter.export2md(nbsdir)
	Documenter.export2readme()
end

# ╔═╡ e2e5d93e-d721-4e71-8e3a-cee7c323ff50
Documenter.showdoc(build)

# ╔═╡ d30b0dc5-8f18-4661-8ecf-147d21480a4d
#hide
data = YAML.load_file("../mkdocs.yml")

# ╔═╡ 9320c981-7d6d-4d26-a1ed-782ec23bda90
#export
mkdocs = OrderedDict("site-name" => "your site name", "nav" => [OrderedDict("Home" => "index.md")],
"theme" => OrderedDict("name" => "material", "palette" => OrderedDict("primary" => "white"), "font" => OrderedDict("text" => "Roboto"), "logo" => "assets/asset image filename", "favicon" => "asset/ asset name", "features" => "navigation.tabs", "markdown_extensions" => ["meta", "toc"], "plugins" => ["search", "tags"]), "extra" => OrderedDict("social" => [OrderedDict("icon" => "fontawsome/brans/twitter", "link" => "twitter link"), OrderedDict("icon" => "fotawsome/brands/linkedin", "link" => "linkedin link")]), "repo_url" => "url to your git repo", "repo_name" => "name of your git repo")

# ╔═╡ e3185d2a-1763-482e-9c5b-b6f9e1aafd84
md"## new"

# ╔═╡ 68e3ea9b-a6a9-4ba6-ac0d-e5c8d4644eae
#export
begin
"""
> new(projname)-> Scaffolds a new project directory with the necessary sub-directories and configuration files.
"""
#TODO: creation of assets, iamges, index.jl is not getting
#created when this function is called from terminal
function new(projname = "newproj")
	subdirs = ["$projname/docs", "nbs", "$projname/src", "$projname/docs/images", "$projname/docs/assets"]
	files = ["$projname/Manifest.toml", "$projname/Project.toml", "$projname/mkdocs.yml", "$projname/nbs/index.jl"]
	map(mkpath, subdirs)
	map(touch, files)
	YAML.write_file(files[3], mkdocs)
end
end

# ╔═╡ 579f18e2-1ada-4f4d-a726-f405763def7a
Documenter.showdoc(new)

# ╔═╡ 02c9cc7c-fd6b-4346-bd36-710d6bc6806d
md"The `new()` function can be used when you want to create a project in a cloned repo. All the files and subdirectories would be blank."

# ╔═╡ 3059e68a-87a3-420f-934c-e2f9bc7499a9
#hide
#new("test1")

# ╔═╡ ae286e66-f1cc-45e7-b271-2a9adf838a0d
testdir = "../test"

# ╔═╡ f481a9f0-d91b-49b6-9bce-48896d6af938
new(joinpath(testdir, "testproj"))

# ╔═╡ 8a115e3d-11d4-4632-8e59-e8c25cb4a295
#noop
@test isdir(joinpath(testdir, "testproj")) == true

# ╔═╡ 785a962e-dfb0-49c0-9679-a6ff0702b001
#noop
@test isdir(joinpath(testdir, "testproj", "nbs")) == true

# ╔═╡ 1de9bc3d-0984-47f3-ad74-adb2290d3889
#noop
@test isdir(joinpath(testdir, "testproj", "docs")) == true

# ╔═╡ 62705a1e-bbca-4bd1-815c-6f00ece58b73
#noop
@test isdir(joinpath(testdir, "testproj", "src")) == true

# ╔═╡ fc17a65b-6223-4f7d-9f75-795781fa81d4
#noop
@test isdir(joinpath(testdir, "testproj", "docs", "images")) == true

# ╔═╡ 5cb8289c-5ae3-4bba-8ac4-0765ee2464a3
#noop
@test isdir(joinpath(testdir, "testproj", "docs", "assets")) == true

# ╔═╡ ccfbfefc-5bc3-45ba-a83d-e9ab2a807311
#noop
@test isfile(joinpath(joinpath(testdir, "testproj"), "Project.toml")) == true

# ╔═╡ d503f5a9-54c9-488b-8ce2-4a81bc0a3cee
#noop
@test isfile(joinpath(joinpath(testdir, "testproj"), "Manifest.toml")) == true

# ╔═╡ 7477e2a6-d2a2-4552-9b03-63daef9c1189
#noop
@test isfile(joinpath(joinpath(testdir, "testproj"), "mkdocs.yml")) == true

# ╔═╡ f29215d9-b898-4540-bf71-6e0c4f192979
#noop
@test isfile(joinpath(joinpath(testdir, "testproj"), "nbs", "index.jl")) == true

# ╔═╡ 89317c3d-5889-4887-893c-aed3bf37b55e
#noop
rm(joinpath(testdir, "testproj"), recursive=true)

# ╔═╡ c96a7fe1-ccba-4499-aa15-95ced972492b
#hide
Export.notebook2script(joinpath("..", "nbs"), joinpath("..", "src"))

# ╔═╡ Cell order:
# ╠═674905ed-75a1-457c-a866-279cdb76d243
# ╠═e3cd4e1f-c7ac-45fd-af9c-1327c0f7a426
# ╠═3e8524c9-1b00-4e30-bb47-524991f3d0f7
# ╠═7a0b47e0-0231-4d3f-b7b8-f6dda83f120d
# ╠═7831e355-bd3b-4333-8200-b8938277af3e
# ╠═2ac37248-81d0-4c8b-99c7-4dd9f400a955
# ╠═2c396254-932f-4431-8f94-08dd2766aa8e
# ╠═17a13720-29ba-11ec-0ba4-83c9ee9ac582
# ╠═1f7d1cb4-7fb9-4f40-9350-70ef40bca38f
# ╠═1c5d5f39-0349-454b-aa22-5feda1287b9d
# ╠═e2e5d93e-d721-4e71-8e3a-cee7c323ff50
# ╠═d30b0dc5-8f18-4661-8ecf-147d21480a4d
# ╠═9320c981-7d6d-4d26-a1ed-782ec23bda90
# ╠═e3185d2a-1763-482e-9c5b-b6f9e1aafd84
# ╠═68e3ea9b-a6a9-4ba6-ac0d-e5c8d4644eae
# ╠═579f18e2-1ada-4f4d-a726-f405763def7a
# ╠═02c9cc7c-fd6b-4346-bd36-710d6bc6806d
# ╠═3059e68a-87a3-420f-934c-e2f9bc7499a9
# ╠═ae286e66-f1cc-45e7-b271-2a9adf838a0d
# ╠═f481a9f0-d91b-49b6-9bce-48896d6af938
# ╠═8a115e3d-11d4-4632-8e59-e8c25cb4a295
# ╠═785a962e-dfb0-49c0-9679-a6ff0702b001
# ╠═1de9bc3d-0984-47f3-ad74-adb2290d3889
# ╠═62705a1e-bbca-4bd1-815c-6f00ece58b73
# ╠═fc17a65b-6223-4f7d-9f75-795781fa81d4
# ╠═5cb8289c-5ae3-4bba-8ac4-0765ee2464a3
# ╠═ccfbfefc-5bc3-45ba-a83d-e9ab2a807311
# ╠═d503f5a9-54c9-488b-8ce2-4a81bc0a3cee
# ╠═7477e2a6-d2a2-4552-9b03-63daef9c1189
# ╠═f29215d9-b898-4540-bf71-6e0c4f192979
# ╠═89317c3d-5889-4887-893c-aed3bf37b55e
# ╠═c96a7fe1-ccba-4499-aa15-95ced972492b


text:
module Utils

function with_log_on_error(f)
    try
        return f()
    catch err
        @error "Error from `$f`" exception = (err, catch_backtrace())
        rethrow()
    end
end

concurrently(functions...; kwargs...) = run_concurrently(collect(functions); kwargs...)

"""
    run_concurrently(thunks; [spawn = false])

Run `thunks` concurrently.
"""
function run_concurrently(functions::AbstractVector; spawn = false)
    tasks = Task[]
    Base.Experimental.@sync begin  # OK for testing
        for f in functions[1:end-1]
            t = if spawn
                Threads.@spawn with_log_on_error(f)
            else
                @async with_log_on_error(f)
            end
            push!(tasks, t)
        end
        global TASKS = tasks  # for debugging
        functions[end]()
    end
end

# This use of `schedule` is wrong but it's handy when debugging implementation
# of concurrent program that is wrong.
function cancel_tasks(tasks = TASKS)
    for t in tasks
        istaskdone(t) && continue
        schedule(t, ErrorException("!!! cancel !!!"); error = true)
    end
    return tasks
end

function random_sleep(spawn::Bool = true)
    if rand() < 0.1
        # no sleep
    elseif spawn && rand(Bool)
        nspins = rand(0:10000)
        for _ in 1:nspins
            GC.safepoint()
            ccall(:jl_cpu_pause, Cvoid, ())
        end
    else
        sleep(rand() / 1_000_000)
    end
end

end  # module


text:
module SimpleDirectMediaLayer

    include("LibSDL2.jl")
    using .LibSDL2

    using ColorTypes

    import Base.unsafe_convert
    export  TTF_Init, TTF_OpenFont, TTF_RenderText_Blended, TTF_SizeText

    mutable struct SDLWindow
        win::Ptr{SDL_Window}
        renderer::Ptr{SDL_Renderer}

        function SDLWindow(w,h,title="SDL Window")
            win = SDL_CreateWindow(title, Int32(100), Int32(100), Int32(w), Int32(h), UInt32(SDL_WINDOW_SHOWN | SDL_WINDOW_INPUT_FOCUS))
            SDL_SetWindowResizable(win,true)

            renderer = SDL_CreateRenderer(win, Int32(-1), UInt32(SDL_RENDERER_ACCELERATED | SDL_RENDERER_PRESENTVSYNC) )

            new(win, renderer)
        end
    end

    function init()
        SDL_GL_SetAttribute(SDL_GL_MULTISAMPLEBUFFERS, 4)
        SDL_GL_SetAttribute(SDL_GL_MULTISAMPLESAMPLES, 4)
        SDL_Init(UInt32(SDL_INIT_VIDEO))
        TTF_Init()
        Mix_OpenAudio(Int32(22050), UInt16(MIX_DEFAULT_FORMAT), Int32(2), Int32(1024) )
    end

    function mouse_position()
        x,y = Int[1], Int[1]
        SDL_GetMouseState(pointer(x), pointer(y))
        x[1],y[1]
    end

    function event()
        ev = SDL_Event(ntuple(i->UInt8(0),56))
        SDL_PollEvent(pointer_from_objref(ev)) == 0 && return nothing

        t = UInt32(0)
        for x in ev._Event[4:-1:1]
            t = t << (sizeof(x)*8)
            t |= x
        end
        evtype = SDL_Event(t)
        evtype == nothing && return nothing

        unsafe_load( Ptr{evtype}(pointer_from_objref(ev)) )
    end

end # module


text:
"""
$(TYPEDEF)
$(TYPEDFIELDS)

"""
struct Float <: Types.AbstractDelayed
    arg::Types.AbstractDelayed
end

Base.float(arg::Types.AbstractDelayed) = Float(arg)


text:
###

const bracketing_error = """The interval [a,b] is not a bracketing interval.
You need f(a) and f(b) to have different signs (f(a) * f(b) < 0).
Consider a different bracket or try fzero(f, c) with an initial guess c.

"""

# pushed this into both __middle and bisection(f,a,b)
@deprecate bisection64(f, a, b) Roots.bisection(f, a, b)

abstract type AbstractBisection <: AbstractBracketing end

"""

    Bisection()

If possible, will use the bisection method over `Float64` values. The
bisection method starts with a bracketing interval `[a,b]` and splits
it into two intervals `[a,c]` and `[c,b]`, If `c` is not a zero, then
one of these two will be a bracketing interval and the process
continues. The computation of `c` is done by `_middle`, which
reinterprets floating point values as unsigned integers and splits
there. This method avoids floating point issues and when the
tolerances are set to zero (the default) guarantees a "best" solution
(one where a zero is found or the bracketing interval is of the type
`[a, nextfloat(a)]`).

When tolerances are given, this algorithm terminates when the midpoint
is approximately equal to an endpoint using absolute tolerance `xatol`
and relative tolerance `xrtol`.

When a zero tolerance is given and the values are not `Float64`
values, this will call the `A42` method.

"""
struct Bisection <: AbstractBisection end  # either solvable or A42
struct BisectionExact <: AbstractBisection end

## tracks for bisection, different, we show bracketing interval
function log_step(l::Tracks, M::AbstractBracketing, state)
    push!(l.xs, state.xn0)
    push!(l.xs, state.xn1) # we store [ai,bi, ai+1, bi+1, ...]
end
function show_tracks(l::Tracks, M::AbstractBracketing)

    xs = l.xs
    n = length(xs)
    for (i,j) in enumerate(1:2:(n-1))
        println(@sprintf("(%s, %s) = (% 18.16f, % 18.16f)", "a_$(i-1)", "b_$(i-1)", xs[j], xs[j+1]))
    end
    println("")
end

## helper function
function adjust_bracket(x0)
    u, v = float.(promote(x0...))
    if u > v
        u,v = v,u
    end
    isinf(u) && (u = nextfloat(u))
    isinf(v) && (v = prevfloat(v))
    u, v
end

# In init_state the bracketing interval is left as (a,b) with
# a,b both finite and of the same sign
function init_state(method::AbstractBisection, fs, x)
    length(x) > 1 || throw(ArgumentError(bracketing_error))

    x0, x1 = adjust_bracket(x) # now finite, right order
    fx0, fx1 = promote(fs(x0), fs(x1))
    sign(fx0) * sign(fx1) > 0 && throw(ArgumentError(bracketing_error))

    _init_state(method, fs, (x0, x1), (fx0, fx1))
end

# assume xs, fxs, have all been checked so that we have a bracket
# gives an  entry for the case where the endpoints are expensive to compute
# as requested in issue #156 by @greimel
function _init_state(method::AbstractBisection, fs, xs, fxs)

    x0, x1 = xs
    fx0, fx1 = fxs

    state = UnivariateZeroState(x1, x0, zero(x1)/zero(x1)*oneunit(x1), [x1],
                                fx1, fx0, fx1, [fx1],
                                0, 2,
                                false, false, false, false,
                                "")

    init_state!(state, method, fs, (x0, x1), (fx0, fx1))
    state
end

function init_state!(state::UnivariateZeroState{T,S}, M::AbstractBisection, fs,
                     xs::Union{Tuple, Vector}) where {T, S}

    x0::T, x1::T = adjust_bracket(xs) # now finite, right order
    fx0::S, fx1::S = promote(fs(x0), fs(x1))
    init_state!(state, M, fs, (x0,x1), (fx0, fx1))

end

function init_state!(state::UnivariateZeroState{T,S}, M::AbstractBisection, fs,
                     xs::Union{Tuple, Vector}, fxs::Union{Tuple, Vector}) where {T, S}
    x0, x1 = xs
    fx0, fx1 = fxs

    sign(fx0) * sign(fx1) > 0 && throw(ArgumentError(bracketing_error))

    if iszero(sign(fx0) * sign(fx1))
        # should be an error -- not bracketing, but we have a zero, so return it.
        if iszero(fx0)
            m, fm = x0, fx0
        else
            m, fm = x1, fx1
        end
        state.f_converged = true
        state.xstar  = m
        state.fxstar =  fm
        #state.xn1 = m
        #state.fxn1 = fm
        return state
    end

    # we need a,b to be same sign, finite
    if sign(x0) * sign(x1) < 0
        m = zero(x1)
        fm::S = fs(m) # iszero(fm) caught in assess_convergence
        incfn(state)
        if sign(fx0) * sign(fm) < 0
            x1, fx1 = m, fm
        elseif  sign(fx0) * sign(fm) > 0
            x0, fx0 = m, fm
        else
            state.message = "Exact zero found"
            state.xstar = m
            state.fxstar = fm
            state.f_converged = true
            state.x_converged = true
            return state
        end
    end

    m = __middle(x0, x1)
    fm = fs(m)
    incfn(state)

    y0, y1, ym = promote(fx0, fx1, fm)

    init_state!(state, x1, x0, [m], y1, y0, [ym])
end

# for Bisection, the defaults are zero tolerances and strict=true
"""
    default_tolerances(M, [T], [S])

For `Bisection` (or `BisectionExact`), when the `x` values are of type `Float64`, `Float32`,
or `Float16`, the default tolerances are zero and there is no limit on
the number of iterations or function evalutions. In this case, the
algorithm is guaranteed to converge to an exact zero, or a point where
the function changes sign at one of the answer's adjacent floating
point values.

For other types, the the `A42` method (with its tolerances) is used.

"""
default_tolerances(M::Union{Bisection, BisectionExact}) = default_tolerances(M,Float64, Float64)
function default_tolerances(::M, ::Type{T}, ::Type{S}) where {M<:Union{Bisection, BisectionExact},T,S}
    xatol = zero(T)
    xrtol = zero(one(T))
    atol = zero(float(one(S))) * oneunit(S)
    rtol = zero(float(one(S))) * one(S)
    maxevals = typemax(Int)
    maxfnevals = typemax(Int)
    strict = true
    (xatol, xrtol, atol, rtol, maxevals, maxfnevals, strict)
end

# find middle of (a,b) with convention that
# * if a, b finite, they are made non-finite
# if a,b of different signs, middle is 0
# middle falls back to a/2 + b/2, but
# for Float64 values, middle is over the
# reinterpreted unsigned integer.
function _middle(x, y)
    a = isinf(x) ? nextfloat(x) : x
    b = isinf(y) ? prevfloat(y) : y
    if sign(a) * sign(b) < 0
        return zero(a)
    else
        __middle(a, b)
    end
end

const FloatNN = Union{Float64, Float32, Float16}

## find middle assuming a,b same sign, finite
## Alternative "mean" definition that operates on the binary representation
## of a float. Using this definition, bisection will never take more than
## 64 steps (over Float64)
__middle(x::Float64, y::Float64) = __middle(Float64, UInt64, x, y)
__middle(x::Float32, y::Float32) = __middle(Float32, UInt32, x, y)
__middle(x::Float16, y::Float16) = __middle(Float16, UInt16, x, y)
## fallback for non FloatNN number types
__middle(x::Number, y::Number) = 0.5*x + 0.5*y

function __middle(T, S, x, y)
    # Use the usual float rules for combining non-finite numbers
    # do division over unsigned integers with bit shift
    xint = reinterpret(S, abs(x))
    yint = reinterpret(S, abs(y))
    mid = (xint + yint) >> 1

    # reinterpret in original floating point
    sign(x+y) * reinterpret(T, mid)
end

## the method converges,
## as we bound between x0, nextfloat(x0) is not measured by eps(), but eps(x0)
function assess_convergence(M::Bisection, state::UnivariateZeroState{T,S}, options) where {T, S}

    assess_convergence(BisectionExact(), state, options) && return true
    x0, x1, xm = state.xn0, state.xn1, first(state.m)

    xtol = max(options.xabstol, max(abs(x0), abs(x1)) * options.xreltol)
    x_converged = x1 - x0 < xtol

    if x_converged
        state.message=""
        state.xstar = xm
        state.x_converged = true
        return true
    end

    fm = first(state.fm)
    ftol = max(options.abstol, max(abs(x0), abs(x1)) * options.reltol)
    f_converged = abs(fm) < ftol

    if f_converged
        state.message = ""
        state.xstar = xm
        state.fxstar = fm
        state.f_converged = f_converged
        return true
    end

    return false
end

# for exact convergence, we can skip some steps
function assess_convergence(M::BisectionExact, state::UnivariateZeroState{T,S}, options) where {T, S}

    state.x_converged && return true

    x0, xm::T, x1 = state.xn0, state.m[1], state.xn1
    y0, ym, y1 = state.fxn0, state.fm[1], state.fxn1

    for (c,fc) in ((x0,y0), (xm,ym), (x1, y1))
        if iszero(fc) || isnan(fc) #|| isinf(fc)
            state.f_converged = true
            state.x_converged = true
            state.xstar = c
            state.fxstar = fc
            return true
        end
    end

    x0 < xm < x1 && return false

    state.xstar = x1
    state.fxstar = ym
    state.x_converged = true
    return true
end

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

function update_state(method::Union{Bisection, BisectionExact}, fs, o::UnivariateZeroState{T,S}, options::UnivariateZeroOptions) where {T<:Number,S<:Number}

    y0 = o.fxn0
    m::T = o.m[1]
    ym::S = o.fm[1] #fs(m)

    if sign(y0) * sign(ym) < 0
        o.xn1, o.fxn1 = m, ym
    else
        o.xn0, o.fxn0 = m, ym
    end

    m  = __middle(o.xn0, o.xn1) # assume a,b have same sign
    fm::S = fs(m)
    o.m[1], o.fm[1] = m, fm
    incfn(o)

    return nothing

end

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

## Bisection has special cases
## for zero tolerance, we have either BisectionExact or A42 methods
## for non-zero tolerances, we have use thegeneral Bisection method
function find_zero(fs, x0, method::M;
                   tracks = NullTracks(),
                   verbose=false,
                   kwargs...) where {M <: Union{Bisection}}

    x = adjust_bracket(x0)
    T = eltype(x[1])
    F = callable_function(fs)
    state = init_state(method, F, x)
    options = init_options(method, state; kwargs...)

    # check if tolerances are exactly 0
    iszero_tol = iszero(options.xabstol) && iszero(options.xreltol) && iszero(options.abstol) && iszero(options.reltol)

    l = (verbose && isa(tracks, NullTracks)) ? Tracks(eltype(state.xn1)[], eltype(state.fxn1)[]) : tracks

    if iszero_tol
        if T <: FloatNN
            return find_zero(F, x, BisectionExact(); tracks=l, verbose=verbose, kwargs...)
        else
            return find_zero(F, x, A42(); tracks=l, verbose=verbose, kwargs...)
        end
    end

    find_zero(method, F, options, state, l)

    verbose && show_trace(method, nothing, state, l)

    state.xstar

end

###################################################
#
## Alefeld, Potra, Shi have two algorithms belosw, one is most efficient, but
## slightly slower than other.
abstract type AbstractAlefeldPotraShi <: AbstractBracketing end

"""
    Roots.A42()

Bracketing method which finds the root of a continuous function within
a provided interval [a, b], without requiring derivatives. It is based
on algorithm 4.2 described in: 1. G. E. Alefeld, F. A. Potra, and
Y. Shi, "Algorithm 748: enclosing zeros of continuous functions," ACM
Trans. Math. Softw. 21, 327–344 (1995), DOI: 10.1145/210089.210111.
Originally by John Travers

"""
struct A42 <: AbstractAlefeldPotraShi end

## put in utils?
@inline isbracket(fa,fb) = sign(fa) * sign(fb) < 0

# f[a, b]
@inline f_ab(a,b,fa,fb) = (fb - fa) / (b-a)

# f[a,b,d]
@inline function f_abd(a,b,d,fa,fb,fd)
    fab, fbd = f_ab(a,b,fa,fb), f_ab(b,d,fb,fd)
    (fbd - fab)/(d-a)
end

# a bit better than a - fa/f_ab
@inline secant_step(a, b, fa, fb) =  a - fa * (b - a) / (fb - fa)

# assume fc != 0
## return a1,b1,d with a < a1 <  < b1 < b, d not there
@inline function bracket(a,b,c, fa, fb, fc)
    if isbracket(fa, fc)
        # switch b,c
        return (a,c,b, fa, fc, fb)
    else
        # switch a,c
        return (c,b,a, fc, fb, fa)
    end
end

# Cubic if possible, if not, quadratic(3)
function take_a42_step(a::T, b, d, ee, fa, fb, fd, fe, k, delta=zero(T)) where {T}

    # if r is NaN or Inf we move on by condition. Faster than checking ahead of time for
    # distinctness
    r = ipzero(a,b,d,ee, fa, fb,fd,fe, delta) # let error and see difference in allcoation?
    (a + 2delta < r < b - 2delta) && return r
    r = newton_quadratic(a,b,d,fa,fb,fd, 3, delta)
end

function ipzero(a::T, b, c, d, fa, fb, fc, fd, delta=zero(T)) where {T}
    Q11 = (c - d)*fc/(fd - fc)
    Q21 = (b - c)*fb/(fc - fb)
    Q31 = (a - b)*fa/(fb - fa)
    D21 = (b - c)*fc/(fc - fb)
    D31 = (a - b)*fb/(fb - fa)
    Q22 = (D21 - Q11)*fb/(fd - fb)
    Q32 = (D31 - Q21)*fa/(fc - fa)
    D32 = (D31 - Q21)*fc/(fc - fa)
    Q33 = (D32 - Q22)*fa/(fd - fa)
    c = a + (Q31 + Q32 + Q33)

    (a + 2delta < c < b - 2delta) && return c

    newton_quadratic(a,b,d,fa,fb,fd, 3, delta)

end

# return c in (a+delta, b-delta)
# adds part of `bracket` from paper with `delta`
function newton_quadratic(a::T, b, d, fa, fb, fd, k::Int, delta=zero(T)) where {T}

    A = f_abd(a,b,d,fa,fb,fd)
    r = isbracket(A,fa) ? b : a

    # use quadratic step; if that fails, use secant step; if that fails, bisection
    if !(isnan(A) || isinf(A)) || !iszero(A)
        B = f_ab(a,b,fa,fb)

        for i in 1:k
            Pr = fa + B * (r-a) +  A * (r-a)*(r-b)
            Prp = (B + A*(2r - a - b))
            r -= Pr / Prp
        end
        if a+2delta < r < b - 2delta
            return r
        end
    end
    # try secant step
    r =  secant_step(a, b, fa, fb)

    if a + 2delta < r < b - 2delta
        return r
    end

    return _middle(a, b) # is in paper r + sgn * 2 * delta

end

# (todo: DRY up?)
function init_state(M::AbstractAlefeldPotraShi, f, xs)
    u, v = promote(float(xs[1]), float(xs[2]))
    if u > v
        u, v = v, u
    end
    fu, fv = promote(f(u), f(v))
    isbracket(fu, fv) || throw(ArgumentError(bracketing_error))
    state = UnivariateZeroState(v, u, zero(v)/zero(v)*oneunit(v), [v, v], ## x1, x0, d, [ee]
                                fv, fu, fv, [fv,fv], ## fx1, fx0, d, [fe]
                                0, 2,
                                false, false, false, false,
                                "")

    init_state!(state, M, f, (u,v), false)
    state
end

# secant step, then bracket for initial setup
function init_state!(state::UnivariateZeroState{T,S}, ::AbstractAlefeldPotraShi, f, xs::Union{Tuple, Vector}, compute_fx=true) where {T, S}

    if !compute_fx
        a, b = state.xn0, state.xn1
        fa, fb = state.fxn0, state.fxn1
    else
        a, b = promote(float(xs[1]), float(xs[2]))
        if a > b
            a, b = b, a
        end
        fa, fb = f(a), f(b)
        state.fnevals = 2
        isbracket(fa, fb) || throw(ArgumentError(bracketing_error))
    end

    c::T = _middle(a, b)
    fc::S = f(c)
    incfn(state)

    a,b,d,fa,fb,fd = bracket(a,b,c,fa,fb,fc)
    ee, fe = d, fd

    init_state!(state, b, a, [d,ee], fb, fa, [fd,fe])
    state.steps = 0
    state.stopped = state.x_converged = state.f_converged = state.convergence_failed = false

    return nothing
end

# for A42, the defaults are reltol=eps(), atol=0; 45 evals and strict=true
# this *basically* follows the tol in the paper (2|u|*rtol + atol)
"""
    default_tolerances(::AbstractAlefeldPotraShi, T, S)

The default tolerances for Alefeld, Potra, and Shi methods are
`xatol=zero(T)`, `xrtol=eps(T)/2`, `atol= zero(S), and rtol=zero(S)`, with
appropriate units; `maxevals=45`, `maxfnevals = Inf`; and `strict=true`.

"""
default_tolerances(M::AbstractAlefeldPotraShi) = default_tolerances(M, Float64, Float64)
function default_tolerances(::AbstractAlefeldPotraShi, ::Type{T}, ::Type{S}) where {T,S}
    xatol = zero(T)
    xrtol = eps(one(T))/2
    atol = zero(float(one(S))) * oneunit(S)
    rtol = zero(float(one(S))) * one(S)
    maxevals = 45
    maxfnevals = typemax(Int)
    strict = true
    (xatol, xrtol, atol, rtol, maxevals, maxfnevals, strict)
end

## convergence is much different here
function check_zero(::AbstractBracketing, state, c, fc)
    if isnan(c)
        state.stopped = true
        #state.xn1 = c
        state.xstar = c
        state.fxstar =fc
        state.message *= "NaN encountered. "
        return true
    elseif isinf(c)
        state.stopped = true
        #state.xn1 = c
        state.xstar = c
        state.fxstar =fc
        state.message *= "Inf encountered. "
        return true
    elseif iszero(fc)
        state.f_converged=true
        state.message *= "Exact zero found. "
        state.xstar = c
        state.fxstar =  fc
#        state.xn1 = c
#        state.fxn1 = fc
        return true
    end
    return false
end

function assess_convergence(method::AbstractAlefeldPotraShi, state::UnivariateZeroState{T,S}, options) where {T,S}

    if state.stopped || state.x_converged || state.f_converged
        ## if isnan(state.xstar)
        ##     state.xstar, state.fxstar = state.xn1, state.fxn1
        ## end
        return true
    end

    if state.steps > options.maxevals
        state.stopped = true
        state.message *= "Too many steps taken. "
        return true
    end

    if state.fnevals > options.maxfnevals
        state.stopped=true
        state.message *= "Too many function evaluations taken. "
        return true
    end

    # check f
    u,fu = choose_smallest(state.xn0, state.xn1, state.fxn0, state.fxn1)
    #u,fu = choose_smallest(u, state.m[1], fu, state.fm[1])

    if abs(fu) <= maximum(promote(options.abstol, abs(u) * oneunit(fu) / oneunit(u) * options.reltol))
        state.f_converged = true
        state.xstar=u
        state.fxstar=fu
        if iszero(fu)

            state.x_converged
            state.message *= "Exact zero found. "
        end
        return true
    end

    a,b = state.xn0, state.xn1
    mx = max(abs(a), abs(b))
    tol = maximum(promote(options.xabstol, mx * options.xreltol, sign(options.xreltol) *   eps(mx)))

    if abs(b-a) <= 2tol
        # use smallest of a,b,m
        state.xstar = u
        state.fxstar = fu
        state.x_converged = true
        return true
    end

    return false
end

## initial step, needs to log a,b,d
function log_step(l::Tracks, M::AbstractAlefeldPotraShi, state, ::Any)
    a, b, c = state.xn0, state.xn1, state.m[1]
    append!(l.xs, extrema((a,b,c)))
    push!(l.xs, a)
    push!(l.xs, b) # we store [ai,bi, ai+1, bi+1, ...] for brackecting methods
end

# Main algorithm for A42 method
function update_state(M::A42, f, state::UnivariateZeroState{T,S}, options::UnivariateZeroOptions) where {T,S}

    a::T,b::T,d::T, ee::T = state.xn0, state.xn1, state.m[1], state.m[2]
    fa::S,fb::S,fd::S,fe::S = state.fxn0, state.fxn1, state.fm[1], state.fm[2]

    an, bn = a, b
    μ, λ = 0.5, 0.7
    tole = max(options.xabstol, max(abs(a),abs(b)) * options.xreltol) # paper uses 2|u|*rtol + atol
    delta = λ * tole

    if state.steps < 1
        c = newton_quadratic(a, b, d, fa, fb, fd, 2)
    else
        c = ipzero(a, b, d, ee, fa, fb, fd, fe)
    end
    fc::S = f(c)
    incfn(state)
    if check_zero(M, state, c, fc)
        return nothing
    end

    ab::T, bb::T, db::T, fab::S, fbb::S, fdb::S = bracket(a,b,c,fa,fb,fc)
    eb::T, feb::S = d, fd

    cb::T = take_a42_step(ab, bb, db, eb, fab, fbb, fdb, feb, delta)
    fcb::S = f(cb)
    incfn(state)
    if check_zero(M, state, cb, fcb)
        # tighten up bracket
        state.xn0, state.xn1, state.m[1]  = ab, bb, db
        state.fxn0, state.fxn1, state.fm[1]= fab, fbb, fdb

        return nothing
    end

    ab,bb,db,fab,fbb,fdb = bracket(ab,bb,cb,fab,fbb,fcb)

    u::T, fu::S = choose_smallest(ab, bb, fab, fbb)
    cb = u - 2 * fu * (bb - ab) / (fbb - fab)
    ch::T = cb
    if abs(cb - u) > 0.5 * (b-a)
        ch = _middle(an, bn)
    end
    fch::S = f(cb)
    incfn(state)
    if check_zero(M, state, ch, fch)
        # tighten up bracket
        state.xn0, state.xn1, state.m[1]  = ab, bb, db
        state.fxn0, state.fxn1, state.fm[1]= fab, fbb, fdb

        return nothing
    end

    ah::T, bh::T, dh::T, fah::S, fbh::S, fdh::S = bracket(ab, bb, ch, fab, fbb, fch)

    if bh - ah < μ * (b - a)
        #a, b, d, fa, fb, fd = ahat, b, dhat, fahat, fb, fdhat # typo in paper
        a, b, d, ee =  ah, bh, dh, db
        fa, fb, fd, fe = fah, fbh, fdh, fdb
    else
        m::T = _middle(ah, bh)
        fm::S = f(m)
        incfn(state)
        ee, fe = dh, fdh
        a, b, d, fa, fb, fd = bracket(ah, bh, m, fah, fbh, fm)
    end
    state.xn0, state.xn1, state.m[1], state.m[2],  = a, b, d, ee
    state.fxn0, state.fxn1, state.fm[1], state.fm[2] = fa, fb, fd, fe

    return nothing
end

####
"""
    Roots.AlefeldPotraShi()

Follows algorithm in "ON ENCLOSING SIMPLE ROOTS OF NONLINEAR
EQUATIONS", by Alefeld, Potra, Shi; DOI:
10.1090/S0025-5718-1993-1192965-2
[link](http://www.ams.org/journals/mcom/1993-61-204/S0025-5718-1993-1192965-2/S0025-5718-1993-1192965-2.pdf). Efficiency
is 1.618. Less efficient, but can be faster than A42() method.

"""
struct AlefeldPotraShi <: AbstractAlefeldPotraShi end

# ## 3, maybe 4, functions calls per step
function update_state(M::AlefeldPotraShi, f, state::UnivariateZeroState{T,S}, options::UnivariateZeroOptions) where {T,S}

    a::T,b::T,d::T = state.xn0, state.xn1, state.m[1]
    fa::S,fb::S,fd::S = state.fxn0, state.fxn1, state.fm[1]

    μ, λ = 0.5, 0.7
    tole = max(options.xabstol, max(abs(a),abs(b)) * options.xreltol) # paper uses 2|u|*rtol + atol
    delta = λ * tole

    c::T = newton_quadratic(a, b, d, fa, fb, fd, 2, delta)
    fc::S = f(c)
    incfn(state)
    check_zero(M, state, c, fc) && return nothing

    a,b,d,fa,fb,fd = bracket(a,b,c,fa,fb,fc)

    c = newton_quadratic(a,b,d,fa,fb,fd, 3, delta)
    fc = f(c)
    incfn(state)
    if check_zero(M, state, c, fc)
        # tighten up bracket
        state.xn0, state.xn1, state.m[1]  = a, b, d
        state.fxn0, state.fxn1, state.fm[1]= fa, fb, fd

        return nothing
    end

    a, b, d, fa, fb, fd = bracket(a, b, c, fa, fb,fc)

    u::T, fu::S = choose_smallest(a, b, fa, fb)
    c = u - 2 * fu * (b - a) / (fb - fa)
    if abs(c - u) > 0.5 * (b - a)
        c = _middle(a, b)
    end
    fc = f(c)
    incfn(state)
    if  check_zero(M, state, c, fc)
        # tighten up bracket
        state.xn0, state.xn1, state.m[1]  = a, b, d
        state.fxn0, state.fxn1, state.fm[1]= fa, fb, fd

        return nothing
    end

    ahat::T, bhat::T, dhat::T, fahat::S, fbhat::S, fdhat::S = bracket(a, b, c, fa, fb, fc)
    if bhat - ahat < μ * (b - a)
        #a, b, d, fa, fb, fd = ahat, b, dhat, fahat, fb, fdhat # typo in paper
        a, b, d, fa, fb, fd = ahat, bhat, dhat, fahat, fbhat, fdhat
    else
        m::T = _middle(ahat, bhat)
        fm::S = f(m)
        incfn(state)
        a, b, d, fa, fb, fd = bracket(ahat, bhat, m, fahat, fbhat, fm)
    end
    state.xn0, state.xn1, state.m[1] = a, b, d
    state.fxn0, state.fxn1, state.fm[1] = fa, fb, fd

    return nothing
end

### Brent
"""
    Roots.Brent()

An implementation of
[Brent's](https://en.wikipedia.org/wiki/Brent%27s_method) (or Brent-Dekker) method.
This method uses a choice of inverse quadratic interpolation or a secant
step, falling back on bisection if necessary.

"""
struct Brent <: AbstractBracketing end

function log_step(l::Tracks, M::Brent, state)
    a,b = state.xn0, state.xn1
    u,v = a < b ? (a,b) : (b,a)
    push!(l.xs, a)
    push!(l.xs, b) # we store [ai,bi, ai+1, bi+1, ...]
end

#
function init_state(M::Brent, f, xs)
    u, v = promote(float(xs[1]), float(xs[2]))
    fu, fv = promote(f(u), f(v))
    isbracket(fu, fv) || throw(ArgumentError(bracketing_error))

    # brent store b as smaller of |fa|, |fb|
    if abs(fu) > abs(fv)
        a, b, fa, fb = u, v, fu, fv
    else
        a, b, fa, fb = v, u, fv, fu
    end

    state = UnivariateZeroState(b, a, zero(b)/zero(b)*oneunit(b), [a, a], ## x1, x0, c, d
                                fb, fa, fb, [fa, one(fa)], ## fx1, fx0, fc, mflag
                                0, 2,
                                false, false, false, false,
                                "")
    state
end

# we store mflag as -1, or +1 in state.m[2]
function init_state!(state::UnivariateZeroState{T,S}, ::Brent, f, xs::Union{Tuple, Vector}) where {T, S}
    u::T, v::T = promote(float(xs[1]), float(xs[2]))
    fu::S, fv::S = promote(f(u), f(v))
    isbracket(fu, fv) || throw(ArgumentError(bracketing_error))

    # brent store b as smaller of |fa|, |fb|
    if abs(fu) > abs(fv)
        a, b, fa, fb = u, v, fu, fv
    else
        a, b, fa, fb = v, u, fv, fu
    end

    init_state!(state, b,a,[a,a], fb,fa, [fa, one(fa)])
    state.steps = 0
    state.stopped = state.x_converged = state.f_converged = state.convergence_failed = false

    return nothing
end

function update_state(M::Brent, f, state::UnivariateZeroState{T,S}, options::UnivariateZeroOptions) where {T,S}
    mflag = state.fm[2] > 0
    a, b, c, d = state.xn0, state.xn1, state.m[1], state.m[2]
    fa, fb, fc = state.fxn0, state.fxn1, state.fm[1]

    # next setp
    s::T = zero(a)
    if  !iszero(fa-fc) && !iszero(fb-fc)
        s =  a * fb * fc / (fa - fb) / (fa - fc) # quad step
        s += b * fa * fc / (fb - fa) / (fb - fc)
        s += c * fa * fb / (fc - fa) / (fc - fb)
        s
    else
        s = secant_step(a,b,fa,fb)
    end
    fs::S = f(s)
    incfn(state)
    check_zero(M, state, s, fs) && return nothing

    # guard step
    u,v = (3a+b)/4, b
    if u > v
        u,v = v, u
    end

    tol = max(options.xabstol, max(abs(b), abs(c), abs(d)) * options.xreltol)
    if !(u < s < v) ||
        (mflag && abs(s - b) >= abs(b-c)/2) ||
        (!mflag && abs(s - b) >= abs(b-c)/2) ||
        (mflag && abs(b-c) <= tol) ||
        (!mflag && abs(c-d) <= tol)
        s = _middle(a, b)
        fs = f(s)
        incfn(state)
        check_zero(M, state, s, fs) && return nothing
        mflag = true
    else
        mflag = false
    end

     d = c
    c,fc = b,fb

    if sign(fa) * sign(fs) < 0
        b, fb = s, fs
    else
        a, fa = s, fs
    end

    if abs(fa) < abs(fb)
        a, b, fa, fb = b, a, fb, fa
    end

    state.xn0, state.xn1, state.m[1], state.m[2] = a, b, c, d
    state.fxn0, state.fxn1, state.fm[1] = fa, fb, fc
    state.fm[2] = mflag ? one(fa) : -one(fa)

    return nothing
end

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

"""

    FalsePosition()

Use the [false
position](https://en.wikipedia.org/wiki/False_position_method) method
to find a zero for the function `f` within the bracketing interval
`[a,b]`.

The false position method is a modified bisection method, where the
midpoint between `[a_k, b_k]` is chosen to be the intersection point
of the secant line with the x axis, and not the average between the
two values.

To speed up convergence for concave functions, this algorithm
implements the 12 reduction factors of Galdino (*A family of regula
falsi root-finding methods*). These are specified by number, as in
`FalsePosition(2)` or by one of three names `FalsePosition(:pegasus)`,
`FalsePosition(:illinois)`, or `FalsePosition(:anderson_bjork)` (the
default). The default choice has generally better performance than the
others, though there are exceptions.

For some problems, the number of function calls can be greater than
for the `bisection64` method, but generally this algorithm will make
fewer function calls.

Examples
```
find_zero(x -> x^5 - x - 1, (-2, 2), FalsePosition())
```
"""
struct FalsePosition{R} <: AbstractBisection end
FalsePosition(x=:anderson_bjork) = FalsePosition{x}()

function update_state(method::FalsePosition, fs, o::UnivariateZeroState{T,S}, options::UnivariateZeroOptions) where {T,S}

    a::T, b::T =  o.xn0, o.xn1
    fa::S, fb::S = o.fxn0, o.fxn1

    lambda = fb / (fb - fa)
    tau = 1e-10                   # some engineering to avoid short moves; still fails on some
    if !(tau < abs(lambda) < 1-tau)
        lambda = 1/2
    end

    x::T = b - lambda * (b-a)
    fx::S = fs(x)
    incfn(o)
    if iszero(fx)
        o.xn1 = x
        o.fxn1 = fx
        o.f_converged = true
        return
    end

    if sign(fx)*sign(fb) < 0
        a, fa = b, fb
    else
        fa = galdino_reduction(method, fa, fb, fx)
    end
    b, fb = x, fx

    o.xn0, o.xn1 = a, b
    o.fxn0, o.fxn1 = fa, fb

    nothing
end

# the 12 reduction factors offered by Galadino
# In RootsTesting.jl, we can see :12 has many more failures.
galdino = Dict{Union{Int,Symbol},Function}(:1 => (fa, fb, fx) -> fa*fb/(fb+fx),
                                           :2 => (fa, fb, fx) -> (fa - fb)/2,
                                           :3 => (fa, fb, fx) -> (fa - fx)/(2 + fx/fb),
                                           :4 => (fa, fb, fx) -> (fa - fx)/(1 + fx/fb)^2,
                                           :5 => (fa, fb, fx) -> (fa -fx)/(1.5 + fx/fb)^2,
                                           :6 => (fa, fb, fx) -> (fa - fx)/(2 + fx/fb)^2,
                                           :7 => (fa, fb, fx) -> (fa + fx)/(2 + fx/fb)^2,
                                           :8 => (fa, fb, fx) -> fa/2,
                                           :9 => (fa, fb, fx) -> fa/(1 + fx/fb)^2,
                                           :10 => (fa, fb, fx) -> (fa-fx)/4,
                                           :11 => (fa, fb, fx) -> fx*fa/(fb+fx),
                                           :12 => (fa, fb, fx) -> (fa * (1-fx/fb > 0 ? 1-fx/fb : 1/2))
)

# give common names
for (nm, i) in [(:pegasus, 1), (:illinois, 8), (:anderson_bjork, 12)]
    galdino[nm] = galdino[i]
end

# from Chris Elrod; https://raw.githubusercontent.com/chriselrod/AsymptoticPosteriors.jl/master/src/false_position.jl
@generated function galdino_reduction(methods::FalsePosition{R}, fa, fb, fx) where {R}
    f = galdino[R]
    quote
        $Expr(:meta, :inline)
        $f(fa, fb, fx)
    end
end

"""
    find_bracket(f, x0, method=A42(); kwargs...)

For bracketing methods returns an approximate root, the last bracketing interval used, and a flag indicating if an exact zero was found as a named tuple.

With the default tolerances, one  of these should be the case: `exact`  is `true` (indicating termination  of the algorithm due to an exact zero  being identified) or the length of `bracket` is less or equal than `2eps(maximum(abs.(bracket)))`. In the `BisectionExact` case, the 2 could  be replaced by 1, as the bracket, `(a,b)` will satisfy  `nextfloat(a) == b `;  the Alefeld,  Potra, and Shi algorithms don't quite have  that promise.

"""
function find_bracket(fs, x0, method::M=A42(); kwargs...) where {M <: Union{AbstractAlefeldPotraShi, BisectionExact}}
    x = adjust_bracket(x0)
    T = eltype(x[1])
    F = callable_function(fs)
    state = init_state(method, F, x)
    options = init_options(method, state; kwargs...)

    # check if tolerances are exactly 0
    iszero_tol = iszero(options.xabstol) && iszero(options.xreltol) && iszero(options.abstol) && iszero(options.reltol)

    find_zero(method, F, options, state, NullTracks())

    (xstar=state.xstar, bracket=(state.xn0, state.xn1), exact=iszero(state.fxstar))

end


text:
function get_NeoHook_loading()
    _C = range(0.0,  0.005, length=3)
    C = [SymmetricTensor{2,3}((x +1.0, x/10, 0.0, 1.0, 0.0, 1.0)) for x in _C]
    return C
end

@testset "NeoHook" begin
    m = NeoHook(μ=76.9, λ=115.0)

    loading = get_NeoHook_loading()
    check_jld2(m, loading, "NeoHook")#, debug_print=false, OVERWRITE_JLD2=true)
    check_tangents_AD(m,loading)
end


text:
import JSON
export readvariable, mapvariable, floats, tones
export createscadparameters, setscadparameter!, writescadparameters
using DataStructures

# create JSON parameter set for scad
function createscadparameters()
  return OrderedDict("fileFormatVersion" => "1", "parameterSets" => OrderedDict())
end

# set a scad parameter
function setscadparameter!(parameters::OrderedDict, setname::String, key::String, value)
  sets = parameters["parameterSets"]
  if !haskey(sets, setname)
    sets[setname] = OrderedDict{String,String}()
  end
  sets[setname][key] = string(value)
end

# write scad parameters to filename
function writescadparameters(parameters::OrderedDict, filename::String)
  open(filename, "w") do io
    write(io, JSON.json(parameters, 2))
  end
end

# read an environment variable
function readvariable(name::String)
  if !haskey(ENV, name)
    return nothing
  end
  return ENV[name]
end

function mapvariable(fn::Function, name::String)
  if !haskey(ENV, name)
    return nothing
  end
  map(fn, split(readvariable(name)))
end

function floats(var)
  return mapvariable(x->parse(Float64, x), var)
end

function tones(var, A)
  return mapvariable(v->tone(v; A=A), var)
end

max_stars_repo_name: isaacwhanson/scad-flute

max_issues_repo_name: isaacwhanson/scad-flute

max_forks_repo_name: isaacwhanson/scad-flute


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

using Markdown
using InteractiveUtils

# ╔═╡ e1a3f480-12e7-11eb-01db-71003ec47f08
md"""
## What *is* an array, actually?

We have seen several examples of arrays.
When you think about an array, you probably think of a chain of boxes with the same shape, stored one after the other, consecutively in memory. And indeed this is the original meaning of the word "array".

However, as a user, we don't actually care about how the data is actually stored! All we really care about is the functionality of an array. Namely, we want the following properties:

- Elements have an integer index (from 1 to `n`);
- An element at position `i` can be accessed using the syntax `A[i]`;
- An element at position `i` can be stored using `A[i] = i` (if the array is mutable)
- We should be able to see the size of an array.
"""

# ╔═╡ 51a660b4-12f3-11eb-2b8a-f5544e650b1c
md"""
Julia provides an **abstraction** of an array called `AbstractArray`. We can define new types of array that "behave like" arrays but are stored differently.
"""

# ╔═╡ 853ac7d0-12f3-11eb-1b1c-993e336cf23a
md"""
## How does indexing work?
"""

# ╔═╡ 94bb1e6c-12f3-11eb-2e18-bb232d2c9cfe
md"""
First let's think about indexing, which is really the key operation on an array.
Let's define a vector:
"""

# ╔═╡ 314ba27a-12ec-11eb-12ce-67650c0020b7
arr = [34,45,67]

# ╔═╡ ab3afe6e-12f3-11eb-196a-79d26b66f681
md"""
We know that we can index into it using square brackets:
"""

# ╔═╡ b32acfda-12f3-11eb-323d-ef91b9ff9f0d
arr[2]

# ╔═╡ b53bc5ba-12f3-11eb-2b09-d38cc266b069
md"""
This is just syntax? What is actually going on? We can see this by looking at the **lowered** form:
"""

# ╔═╡ 3a202684-12ec-11eb-26fc-df7b7484e671
Meta.@lower arr[2]

# ╔═╡ c225d324-12f3-11eb-23d7-45548e7b32b6
md"""
We see that `arr[2]` is just "syntactic sugar" for calling the `getindex` function as `getindex(arr, 2)`; in other words, it's just a nice short and suggestive way to write the same thing.
"""

# ╔═╡ d92c5aa2-12f3-11eb-286e-8769ecdf03cb
md"""
So to implement our own objects that behave like arrays, we will define a new type and extend `getindex` on that type.
"""

# ╔═╡ f10cf87a-12f3-11eb-1104-0186b9e43757
md"""
## Example: A zero-based vector
"""

# ╔═╡ f7cadd08-12f3-11eb-1882-8dd4ce8ee2c2
md"""
There are interminable discussions online about whether zero-based indexing (i.e. array indices starting at 0) or one-based indexing (starting at 1) is superior. The answer is that *neither* is -- each is more suitable in different contexts.

Fortunately Julia is flexible enough to allow you to use *any* type of indexing that suits you -- see e.g. the `OffsetArrays.jl` package.

Let's see how to implement the simplest case, namely a zero-based vector.

So we define a new type and extend `getindex` on it. In order to use the full Julia machinery we need to extend a couple of other functions too:
"""

# ╔═╡ 1ac0095e-12eb-11eb-3df2-e98d0ff66dd5
struct ZeroFirstVector{T} <: AbstractVector{T}
	array::Vector{T}
end

# ╔═╡ a66afe14-12eb-11eb-3afc-fdac2325ff84
Base.size(vec::ZeroFirstVector) = size(vec.array)

# ╔═╡ e86c53fa-12eb-11eb-0ee4-41455a218d84
Base.axes(vec::ZeroFirstVector) = (0:length(vec.array)-1,)

# ╔═╡ 69b10842-12eb-11eb-381a-0bda698a59ed
Base.getindex(vec::ZeroFirstVector, i::Int) = vec.array[i+1]

# ╔═╡ 9c57e4aa-12eb-11eb-16fb-c7a18c238192
zarr = ZeroFirstVector([1, 3, 4])

# ╔═╡ 3cf51826-12f4-11eb-2405-db9616fa3152
md"""
Now we are able to use zero-based indexing!:
"""

# ╔═╡ fed5e2e4-12eb-11eb-169f-8b050f9eea4a
zarr[0]

# ╔═╡ 464c22f2-12f4-11eb-0b3c-872621c385f6
md"""
## Modifying elements: `setindex!`
"""

# ╔═╡ 4dec992e-12f4-11eb-1f20-b96632769bed
md"""
In a similar way to `getindex`, modifying elements of an array uses the mutating function `setindex!`:
"""

# ╔═╡ 5b8f3508-12ec-11eb-0a31-61ddf796ac74
Meta.@lower arr[2] = 3

# ╔═╡ 0c8ebad2-12ec-11eb-03f1-c5b9c79bb4da
Base.setindex!(vec::ZeroFirstVector, val, i::Int) = (vec.array[i+1] = val)

# ╔═╡ 8438424c-12ec-11eb-1de2-0bb3924ec78c
zarr[0] = 55

# ╔═╡ 9cee3206-12ec-11eb-09c7-1f0ebad11b8e
zarr

# ╔═╡ b0ab2844-12ec-11eb-1908-0318510b3292
md"""
### Julia gives you added bonuses here:
"""

# ╔═╡ bc8154b8-12ec-11eb-1fcc-4765fcd0796d
zarr[0:1]

# ╔═╡ 8c0c8188-12f4-11eb-1576-23958903a69a
@which zarr[0:1]

# ╔═╡ 7425ef6e-12f4-11eb-21c3-19c221c1ebdb
md"""
We never defined how to index using a *range*, only an integer. Julia has a generic fallback mechanism that understands how to do this once you have implemented `getindex` for a single element.

Similarly:
"""

# ╔═╡ c67fbbbc-12ec-11eb-32be-2ff720964fe6
zarr[ [2, 1, 0] ]

# ╔═╡ d13a65b6-12ec-11eb-2996-37987cfae782
zarr[0:1] .= [7,8];

# ╔═╡ dbb95f42-12ec-11eb-01d5-df2e21affd86
zarr

# ╔═╡ e6428aa6-12ec-11eb-06d0-95404a8a4039
zarr[ [2, 1, 0] ] .= [1, 2, 3];

# ╔═╡ e21658ec-12ef-11eb-06fd-b1d282edf7b6
zarr

# ╔═╡ b405b9ac-12f4-11eb-12fa-f7575ff1fa73
md"""
Various other functions also automatically work, e.g.:
"""

# ╔═╡ e57329a2-12ef-11eb-165d-51f080c73652
sort!(zarr)

# ╔═╡ e8b34188-12ef-11eb-2d7f-fde100b582d1
md"""
### Some more examples

"""

# ╔═╡ bd6bfd10-12f4-11eb-17fe-6f9ce83f0ce8
md"""
We have already seen some other examples of array types in the course, e.g. ranges:
"""

# ╔═╡ f0d952da-12ef-11eb-2452-471b184bc529
r = 53:83

# ╔═╡ 0fc7e3e6-12f0-11eb-1b49-e37196154ec6
size(r)

# ╔═╡ 03b23fe0-12f0-11eb-0970-33b3b973e902
r[20]

# ╔═╡ 00088cda-12f0-11eb-36b7-9b900fa6d0bf
typeof(r)

# ╔═╡ f6024924-12ef-11eb-1d6c-3ffff2546914
supertypes(typeof(r))

# ╔═╡ cdf3e532-12f4-11eb-1ce0-873bcf0dea58
md"""
Let's also look at what happens if we reshape an array:
"""

# ╔═╡ f9bf5200-12ef-11eb-30ee-9b3cbbccc886
A = reshape(1:12, 3, 4)

# ╔═╡ 07c1cdb6-12f1-11eb-0b31-15b5c88caac6
typeof(A)

# ╔═╡ 09d69238-12f1-11eb-1d9f-f72e8c6c14d9
supertypes(typeof(A))

# ╔═╡ ff598e1c-12f0-11eb-24f6-bdaa74db2ad7
size(A)

# ╔═╡ 2f23dd4e-12f0-11eb-01fc-53f66d2a55e4
Aᵀ = transpose(A)

# ╔═╡ ea8e0794-12f0-11eb-2194-bdc46224dc69
size(Aᵀ)

# ╔═╡ 2a76a87a-12f1-11eb-3e58-e5f919ba81e8
@which Aᵀ[2,3]

# ╔═╡ 892e7c40-12f0-11eb-2871-0d308063f025
md"The `getindex` method for `Transpose` is defined [here](https://github.com/JuliaLang/julia/blob/master/stdlib/LinearAlgebra/src/adjtrans.jl#L190)"

# ╔═╡ da41895c-12f4-11eb-2ee4-f1ffc4344f61
md"""
## Conclusions
"""

# ╔═╡ dfac92ea-12f4-11eb-1253-b74be632de64
md"""
Summarising, we see that we -- and Julia -- can think of an array as "anything that behaves like an array", i.e. which you can index into.

What does it mean for something to behave like an array? In the end it means that we *could* "materialise" it to produce an actual array that did live as data stored linearly in memory, and our program would behave in exactly the same way.

But Julia provides a mechanism by which we can store less data, or store it in a more convenient way, and have indexing be almost as efficient.
"""

# ╔═╡ Cell order:
# ╟─e1a3f480-12e7-11eb-01db-71003ec47f08
# ╟─51a660b4-12f3-11eb-2b8a-f5544e650b1c
# ╟─853ac7d0-12f3-11eb-1b1c-993e336cf23a
# ╟─94bb1e6c-12f3-11eb-2e18-bb232d2c9cfe
# ╠═314ba27a-12ec-11eb-12ce-67650c0020b7
# ╟─ab3afe6e-12f3-11eb-196a-79d26b66f681
# ╠═b32acfda-12f3-11eb-323d-ef91b9ff9f0d
# ╟─b53bc5ba-12f3-11eb-2b09-d38cc266b069
# ╠═3a202684-12ec-11eb-26fc-df7b7484e671
# ╟─c225d324-12f3-11eb-23d7-45548e7b32b6
# ╟─d92c5aa2-12f3-11eb-286e-8769ecdf03cb
# ╟─f10cf87a-12f3-11eb-1104-0186b9e43757
# ╟─f7cadd08-12f3-11eb-1882-8dd4ce8ee2c2
# ╠═1ac0095e-12eb-11eb-3df2-e98d0ff66dd5
# ╠═a66afe14-12eb-11eb-3afc-fdac2325ff84
# ╠═e86c53fa-12eb-11eb-0ee4-41455a218d84
# ╠═69b10842-12eb-11eb-381a-0bda698a59ed
# ╠═9c57e4aa-12eb-11eb-16fb-c7a18c238192
# ╟─3cf51826-12f4-11eb-2405-db9616fa3152
# ╠═fed5e2e4-12eb-11eb-169f-8b050f9eea4a
# ╟─464c22f2-12f4-11eb-0b3c-872621c385f6
# ╟─4dec992e-12f4-11eb-1f20-b96632769bed
# ╠═5b8f3508-12ec-11eb-0a31-61ddf796ac74
# ╠═0c8ebad2-12ec-11eb-03f1-c5b9c79bb4da
# ╠═8438424c-12ec-11eb-1de2-0bb3924ec78c
# ╠═9cee3206-12ec-11eb-09c7-1f0ebad11b8e
# ╟─b0ab2844-12ec-11eb-1908-0318510b3292
# ╠═bc8154b8-12ec-11eb-1fcc-4765fcd0796d
# ╠═8c0c8188-12f4-11eb-1576-23958903a69a
# ╟─7425ef6e-12f4-11eb-21c3-19c221c1ebdb
# ╠═c67fbbbc-12ec-11eb-32be-2ff720964fe6
# ╠═d13a65b6-12ec-11eb-2996-37987cfae782
# ╠═dbb95f42-12ec-11eb-01d5-df2e21affd86
# ╠═e6428aa6-12ec-11eb-06d0-95404a8a4039
# ╠═e21658ec-12ef-11eb-06fd-b1d282edf7b6
# ╟─b405b9ac-12f4-11eb-12fa-f7575ff1fa73
# ╠═e57329a2-12ef-11eb-165d-51f080c73652
# ╟─e8b34188-12ef-11eb-2d7f-fde100b582d1
# ╟─bd6bfd10-12f4-11eb-17fe-6f9ce83f0ce8
# ╠═f0d952da-12ef-11eb-2452-471b184bc529
# ╠═0fc7e3e6-12f0-11eb-1b49-e37196154ec6
# ╠═03b23fe0-12f0-11eb-0970-33b3b973e902
# ╠═00088cda-12f0-11eb-36b7-9b900fa6d0bf
# ╠═f6024924-12ef-11eb-1d6c-3ffff2546914
# ╟─cdf3e532-12f4-11eb-1ce0-873bcf0dea58
# ╠═f9bf5200-12ef-11eb-30ee-9b3cbbccc886
# ╠═07c1cdb6-12f1-11eb-0b31-15b5c88caac6
# ╠═09d69238-12f1-11eb-1d9f-f72e8c6c14d9
# ╠═ff598e1c-12f0-11eb-24f6-bdaa74db2ad7
# ╠═2f23dd4e-12f0-11eb-01fc-53f66d2a55e4
# ╠═ea8e0794-12f0-11eb-2194-bdc46224dc69
# ╠═2a76a87a-12f1-11eb-3e58-e5f919ba81e8
# ╟─892e7c40-12f0-11eb-2871-0d308063f025
# ╟─da41895c-12f4-11eb-2ee4-f1ffc4344f61
# ╟─dfac92ea-12f4-11eb-1253-b74be632de64

max_stars_repo_name: mathijsvdv/ComputationalThinking

max_issues_repo_name: mathijsvdv/ComputationalThinking

max_forks_repo_name: mathijsvdv/ComputationalThinking