text:
"""
    Given a Mesh, a BoundedFluidQueue, a grid spacing hx, and order construct
a DiscretisedFluidQueue as specified by Mesh, and a Generator. By default
also constructs an initial condition which approximates a point mass at x=eps(),
in phase i=1
"""
function build_discretised_model(type::Type{<:Mesh{T}},model,hx,order) where T
    nodes = range(0.0,model.b; length=Int(round(model.b/hx))+1)
    mesh = type(nodes,order)
    dq = DiscretisedFluidQueue(model,mesh)
    B = build_full_generator(dq)
    d0 = interior_point_mass(eps(),1,dq)
    return d0, B, dq
end

t = StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}}
# some functions to help buid approximations en masse
dg_model(order) = build_discretised_model(DGMesh{t},model,h,order)
order1_model(order) = build_discretised_model(DGMesh{t},model,h/order,1)
qbdrap_model(order) = build_discretised_model(FRAPMesh{t},model,h,order)

"""
Build and save all of the approximations as specified by `orders` and `approx_types`.
`orders` is a vector of odd integers. `approx_types` is a vector of symbols. If `:dg`
appears twice in `approx_types` then a flux limiter is applied to the second one.

`d0_map` is a function which take a `DiscretisedFluidQueue` and maps it to an initial
condition.

`a_map` is a mapping to apply to the approximation which returns the information of interest.
For example, `map_to_t(d0,B,approx_type,t,h=0.005,limit=false,kwargs...)`.

Construct a Dict of `models` with keys `orders`. `models[o]` is another Dict
with keys `approx_types` and values as follows:
    "pdf"=>the_pdf,
    "cdf"=>the_cdf,
    "coeffs_mapped"=>dt,
    "coeffs_0"=>d0,
    "dq"=>dq,
    "B"=>B,
"""
function make_approximations(orders,approx_types,d0_map,a_map=(d0...)->d0[1])
    models = Dict{Any,Dict{Any,Dict}}()
    for (c_o,o) in enumerate(orders)
        tmp_dict = Dict()
        dg_counter = 0
        for (c_m,m) in enumerate(approx_types)

            ~,B,dq = @eval $(Symbol(m,"_model"))($o)

            d0 = d0_map(dq)

            (string(m)=="dg")&&(dg_counter+=1)
            # if two dg methods are specified, use a limiter with the second one
            (dg_counter==2) ? (l=true) : (l=false)

            dt = a_map(d0,B,m,transient_time,0.005,l)

            tmp_dict[string(m,c_m)] = Dict(
                "coeffs_mapped"=>dt,
                "coeffs_0"=>d0,
                "dq"=>dq,
                "B"=>B.B,
            )
        end
        models[o] = tmp_dict
    end
    return models
end

"""
Automatically save the output of `evalate!` under the function name (as a string) as the key.
"""
macro evaluate!(models,funs,grid,_phases)
    expr = quote
        evaluate!(
            $(esc(models)),
            $(esc(funs)),
            $(esc(grid)),
            $(esc(_phases)),
            $(string.(esc(funs).args[1].args)),
        )
    end
    return expr
end

"""
Evaluate the functions in `funs` for all of the `models` at the x-values in `grid`
and phases in `_phases`.

If `f` is in `funs` then we evaluate `f(models[o][m]).(grid[:],_phases[:]')` where
`m` and `o` are keys specifying the approximation type and order, respectively.
The output is stored in `models[o][m][fun_name]`.
"""
function evaluate!(models,funs::Tuple,grid::AbstractVector,
    _phases::AbstractVector{Int},fun_names::Union{AbstractVector{String},Tuple})

    orders = sort([l for l in keys(models)])
    k = sort([k for k in keys(first(models)[2])])

    # all the models are the same, so just grab one of them
    model = models[orders[1]][k[1]]["coeffs_mapped"].dq.model

    for (c_o,o) in enumerate(orders)
        for (c_m,m) in enumerate(k)
            for (c_f,f) in enumerate(funs)
                temp_fun = f(models[o][m]["coeffs_mapped"])
                models[o][m][fun_names[c_f]*"_evaluated"] = temp_fun.(grid[:],_phases[:]')
                models[o][m]["boundary_masses"] = models[o][m]["coeffs_mapped"][[1:N₋(model);end-N₊(model)+1:end]]
            end
        end
    end
    return nothing
end

"""
Evaluate `evaluate(models,funs,grid,phases)` where `phases` are all the phases in the model
"""
macro evaluate!(models,funs,grid)
    expr = quote
        evaluate!(
            $(esc(models)),
            $(esc(funs)),
            $(esc(grid)),
            $(string.(esc(funs).args[1].args)),
        )
    end
    return expr
end
"""
Evaluate `evaluate(models,funs,grid,phases)` where `phases` are all the phases in the model
"""
function evaluate!(models,funs::Tuple,grid::AbstractVector,fun_names::Union{AbstractVector{String},Tuple})
    o = sort([l for l in keys(models)])
    k = keys(models[o[1]])
    model = models[o[1]][k[1]]["coeffs_mapped"].dq.model
    evaluate!(models,funs,grid,1:n_phases(model),fun_names)
    return nothing
end
"""
Call evaluate(sims::Simulation,funs::Tuple,grid::AbstractVector,_phases::AbstractVector{Int},fun_names::Union{AbstractVector{String},Tuple})
where `fun_names` are the names of the functions supplied in `funs`
"""
macro evaluate(sims,funs,grid,_phases)
    # esc(funs).args[1].args -- esc() returns an expression, the first .args gets the args of the output of esc()
    # the second .args return the function names in `funs` as a vector of symbols
    return quote
        evaluate($(esc(sims)),$(esc(funs)),$(esc(grid)),$(esc(_phases)),$(string.(esc(funs).args[1].args)))
    end
end
"""
Evaluate the functions in `funs` for the `sim` at the x-values in `grid`
and phases in `_phases`.

If `f` is in `funs` then we evaluate `f(models[o][m]).(grid[:],_phases[:]')` where
`m` and `o` are keys specifying the approximation type and order, respectively.
The output is stored in `models[o][m][fun_name]`.
"""
function evaluate(sims::Simulation,funs::Tuple,grid::AbstractVector,_phases::AbstractVector{Int},fun_names::Union{AbstractVector{String},Tuple})
    evaluated = Dict{String,Any}()
    for (c_f,f) in enumerate(funs)
        temp_fun = f(sims)
        evaluated[fun_names[c_f]*"_evaluated"] = temp_fun.(grid[:],_phases[:]')
        evaluated["boundary_masses"] = sim_point_masses(sims)
    end
    return evaluated
end
macro evaluate(sims,funs,grid)
    return quote
        evaluate($(esc(sims)),$(esc(funs)),$(esc(grid)),$(string.(esc(funs).args[1].args)))
    end
end
evaluate(sims::Simulation,funs::Tuple,grid::AbstractVector,fun_names::Union{AbstractVector{String},Tuple}) =
    evaluate(sims,funs,grid,1:n_phases(sims.model),fun_names)

"""
evaluate `cell_probs(model).(centres[:],_phases[:]')` for all models as specified in the Dict
`models` (as outpud from `make_approximations()`)
"""
function evaluate_cell_probs!(models,nodes,_phases)
    dict_names = sort([k for k in keys(first(models)[2])])
    orders = sort([l for l in keys(models)])
    centres = (nodes[1:end-1]+nodes[2:end])./2
    for (c_o,o) in enumerate(orders)
        for (c_m,m) in enumerate(dict_names)
            if m[1:end-1]!="order1"
                models[o][m]["cell_probs_evaluated"] = cell_probs(models[o][m]["coeffs_mapped"]).(centres[:],_phases[:]')
            else
                # need to map the coefficients of the order1 method to have the same cell_width
                # as the dg method.
                d = SFMDistribution(
                    permute_order1_coeffs(models[o][m]["coeffs_mapped"].coeffs,
                        models[o][m]["coeffs_mapped"].dq.model,
                        o,
                    ),
                    models[o]["dg1"]["coeffs_mapped"].dq
                )
                models[o][m]["cell_probs_evaluated"] = cell_probs(d).(centres[:],_phases[:]')
            end
        end
    end
    return nothing
end

"""
Simulate, `n_sim` times the `BoundedFluidQueue` as specified by `model` until the `stopping_time` with
initial points specified by `init` where rows of init are samples of [X0, φ0]
"""
function simulate_model(model,stopping_time,n_sims,init,pth,rng::Random.AbstractRNG=Randon.default_rng())
    SFM0 = init(n_sims,rng)
    sims = simulate(model,stopping_time,SFM0,rng)
    df = DataFrame(t=sims.t,φ=sims.φ,X=sims.X,n=sims.n)
    CSV.write(pth*"/data/sim.csv",df)
    return sims
end

"""
compute the values of the point masses as approximated by the data in sims.
"""
function sim_point_masses(sims::Simulation)
    pm_evaluated = zeros(N₋(sims.model)+N₊(sims.model))
    c = 0
    n_sim = length(sims.t)
    for i in 1:n_phases(sims.model)
        if DiscretisedFluidQueues._has_left_boundary(sims.model.S,i)
            c += 1
            pm_evaluated[c] = sum(sims.X[sims.φ.==i]==0.0)/n_sim
        end
    end
    c = 0
    for i in 1:n_phases(sims.model)
        if DiscretisedFluidQueues._has_right_boundary(model.S,i)
            c += 1
            pm_evaluated[c] = sum(sims.X[sims.φ.==i]==sims.model.b)/n_sim
        end
    end
    return pm_evaluated
end

function bootstrap_error_ci(f1_evaluated,f1_pm,boot_samples_evaluated,boot_samples_pm_evaluated,grid,error_metric,p)
    boot_errors = Vector{Float64}(undef,size(boot_samples_evaluated,1))
    for n in 1:size(boot_samples_evaluated,1)
        boot_errors[n] = error_metric(f1_evaluated,f1_pm,boot_samples_evaluated[n,:,:],boot_samples_pm_evaluated[n,:,:],grid)
    end
    return (quantile!(boot_errors,p)..., boot_errors)
end

function analysis(err_funs,models,truth,grid,p,titles,pth; w=true)
    dict_names = sort([k for k in keys(first(models)[2])])
    df_names = copy(dict_names)
    if (sum(df_names.=="dg")>1)
        df_names[findlast(df_names.=="dg")] = "dg_limiter"
    end
    err_names = string.(err_funs[i][1] for i in 1:length(err_funs))
    orders = sort([k for k in keys(models)])
    for (c_f,f) in enumerate(err_funs)
        error_metric = f[1]
        apply_to = f[2]
        err_data = Matrix{Float64}(undef,length(orders),length(dict_names))
        errs_ci_lwr = Matrix{Float64}(undef,length(orders),length(dict_names))
        errs_ci_upr = Matrix{Float64}(undef,length(orders),length(dict_names))
        (w)&&(specific_path=mkpath(pth*"/data/bootstrap_errors/"*err_names[c_f]*"/"*apply_to*"/"))
        for (c_o,o) in enumerate(orders)
            for (c_m,m) in enumerate(dict_names)
                err_data[c_o,c_m] = error_metric(models[o][m][apply_to*"_evaluated"],models[o][m]["boundary_masses"],truth[apply_to*"_evaluated"],truth["boundary_masses"],grid)
                errs_ci_lwr[c_o,c_m], errs_ci_upr[c_o,c_m], bootstrap_data =
                    bootstrap_error_ci(
                        models[o][m][apply_to*"_evaluated"],models[o][m]["boundary_masses"],
                        truth["bootstrap_samples"][apply_to],truth["bootstrap_samples"]["boundary_masses"],
                        grid,error_metric,p,
                    )#error_metric(models[o][m][apply_to*"_evaluated"],models[o][m]["boundary_masses"],truth[apply_to*"_ci_lwr"],truth["boundary_masses"*"_ci_lwr"],grid)
                # errs_ci_upr[c_o,c_m]  = error_metric(models[o][m][apply_to*"_evaluated"],models[o][m]["boundary_masses"],truth[apply_to*"_ci_upr"],truth["boundary_masses"*"_ci_upr"],grid)
                (w)&&JSON.write(specific_path*"/order_"*string(o)*"_model_"*m*".json",JSON.json(bootstrap_data))
            end
        end
        err_data = DataFrame(err_data,df_names)
        write_errs(err_data,pth,err_names[c_f])
        errs_ci_lwr = DataFrame(errs_ci_lwr,df_names)
        write_errs(errs_ci_lwr,pth,err_names[c_f]*"_ci_lwr")
        errs_ci_upr = DataFrame(errs_ci_upr,df_names)
        write_errs(errs_ci_upr,pth,err_names[c_f]*"ci_upr")
        (w)&&make_err_plot(models,err_data,errs_ci_lwr,errs_ci_upr,pth*"/figs/"*err_names[c_f],titles[c_f])
    end
    return nothing
end

function make_err_plot(models,err_data,errs_ci_lwr,errs_ci_upr,pth,title)
    orders = sort([l for l in keys(models)])
    q = plot()
    for (c_m,m) in enumerate(names(err_data))
        lwr = err_data[:,m]-errs_ci_lwr[:,m]
        upr = errs_ci_upr[:,m]-err_data[:,m]
        plot!(q,orders,err_data[:,m],ribbon=(lwr,upr),label=m,color=c_m)
    end
    plot!(q,xlabel="Order",ylabel="log₁₀ Error",title=title,legend=:outertopright)
    savefig(q,pth*".svg")
    return nothing
end

# function compute_errors(orders,models,truth,n_evals)
#     L1_cdf_errors = DataFrame(DG = Vector{Union{Missing, Float64}}(undef,0), DG_Limit = Vector{Union{Missing, Float64}}(undef,0), Order_1 = Vector{Union{Missing, Float64}}(undef,0), QBDRAP = Vector{Union{Missing, Float64}}(undef,0))
#     ks_errors = DataFrame(DG = Vector{Union{Missing, Float64}}(undef,0), DG_Limit = Vector{Union{Missing, Float64}}(undef,0), Order_1 = Vector{Union{Missing, Float64}}(undef,0), QBDRAP = Vector{Union{Missing, Float64}}(undef,0))
#     L1_pdf_errors = DataFrame(DG = Vector{Union{Missing, Float64}}(undef,0), DG_Limit = Vector{Union{Missing, Float64}}(undef,0), Order_1 = Vector{Union{Missing, Float64}}(undef,0), QBDRAP = Vector{Union{Missing, Float64}}(undef,0))
#     L2_pdf_errors = DataFrame(DG = Vector{Union{Missing, Float64}}(undef,0), DG_Limit = Vector{Union{Missing, Float64}}(undef,0), Order_1 = Vector{Union{Missing, Float64}}(undef,0), QBDRAP = Vector{Union{Missing, Float64}}(undef,0))

#     k = keys(models[minimum(keys(models))])

#     for (c_o,o) in enumerate(orders)
#         L1_cdf_error_row = Vector{Union{Missing, Float64}}(undef,length(k))
#         ks_error_row = Vector{Union{Missing, Float64}}(undef,length(k))
#         L1_pdf_error_row = Vector{Union{Missing, Float64}}(undef,length(k))
#         L2_pdf_error_row = Vector{Union{Missing, Float64}}(undef,length(k))
#         for (c_m,m) in enumerate(sort([i for i in k]))
#             L1_cdf_error_row[c_m] = log10(DiscretisedFluidQueues.Lp(x->models[o][m]["cdf"](x,1),x->truth.cdf(x,1),range(-eps(),10+eps(),length=n_evals)) +
#                                                             DiscretisedFluidQueues.Lp(x->models[o][m]["cdf"](x,2),x->truth.cdf(x,2),range(-eps(),10+eps(),length=n_evals)))
#             ks_error_row[c_m] = log10(DiscretisedFluidQueues.kolmogorov_smirnov(x->models[o][m]["cdf"](x,1),x->truth.cdf(x,1),range(-eps(),10+eps(),length=n_evals)) +
#                                                             DiscretisedFluidQueues.kolmogorov_smirnov(x->models[o][m]["cdf"](x,2),x->truth.cdf(x,2),range(-eps(),10+eps(),length=n_evals)))
#             L1_pdf_error_row[c_m] = log10(DiscretisedFluidQueues.Lp(x->models[o][m]["pdf"](x,1),x->truth.pdf(x,1),range(-eps(),10+eps(),length=n_evals)) +
#                                                             DiscretisedFluidQueues.Lp(x->models[o][m]["pdf"](x,2),x->truth.pdf(x,2),range(-eps(),10+eps(),length=n_evals)))
#             L2_pdf_error_row[c_m] = log10(DiscretisedFluidQueues.Lp(x->models[o][m]["pdf"](x,1),x->truth.pdf(x,1),range(-eps(),10+eps(),length=n_evals),2) +
#                                                             DiscretisedFluidQueues.Lp(x->models[o][m]["pdf"](x,2),x->truth.pdf(x,2),range(-eps(),10+eps(),length=n_evals),2))
#         end
#         push!(L1_cdf_errors,L1_cdf_error_row)
#         push!(ks_errors,ks_error_row)
#         push!(L1_pdf_errors,L1_pdf_error_row)
#         push!(L2_pdf_errors,L2_pdf_error_row)
#     end
#     for e in (L1_cdf_errors,ks_errors,L1_pdf_errors,L2_pdf_errors)
#         for i in 1:size(e,1), j in size(e,2)
#             if ismissing(e[i,j])||isnan(e[i,j])
#                 e[i,j]=missing
#             end
#         end
#     end
#     return (L1_cdf_errors,ks_errors,L1_pdf_errors,L2_pdf_errors)
# end

function write_errs(err,pth,err_names)
    file = pth*"/data/"*string(err_names)
    CSV.write(file*".csv",err)
    return nothing
end

std_plot!(p,fun,args...; kwargs...) = plot!(
        p,
        fun,
        _xlims[1], _xlims[2];
        # ylim=(-0.1,2.1),
        xticks=false,
        yticks=false,
        grid=false,
        label=false,
        kwargs...,
    )

function make_df_plots!(models,truth,map_name,grid,_phases,pth,_xlims,_ylims,_xticks,_yticks)
    orders = sort([l for l in keys(models)])
    k = keys(models[orders[1]])
    p1 = plot(layout = (length(k),length(orders)))
    t = keys(truth)
    for (c_o,o) in enumerate(orders)
        for (c_m,m) in enumerate(sort([i for i in k]))
            for (c_i,i) in enumerate(_phases)
                (map_name*"_evaluated"∈t)&&(y_vals = truth[map_name*"_evaluated"][:,c_i])
                if (map_name*"_ci_lwr")∈t
                    lwr = y_vals-truth[map_name*"_ci_lwr"][:,c_i]
                    upr = truth[map_name*"_ci_upr"][:,c_i]-y_vals
                    plot!(p1.layout.grid[c_m,c_o],grid,y_vals,ribbon=(lwr,upr),color=i+maximum(_phases),label=false,xticks=false,alpha=0.5)
                elseif (map_name∈t)
                    plot!(p1.layout.grid[c_m,c_o],grid,y_vals,color=i+maximum(_phases),label=false,xticks=false,alpha=0.5)
                end
                y_vals = models[o][m][map_name*"_evaluated"][:,c_i]
                plot!(p1.layout.grid[c_m,c_o],grid,y_vals,color=i,xticks=false,label=false,alpha=0.5)
            end
            if c_o==1
                plot!(p1.layout.grid[c_m,c_o],ylabel=uppercase(string(m)[1:end-1]))
                plot!(p1.layout.grid[c_m,c_o],yticks=_yticks)
            end
            if c_m==1
                plot!(p1.layout.grid[c_m,c_o],title=string(o))
            end
            if c_m==length(k)
                plot!(p1.layout.grid[c_m,c_o],xticks=_xticks)
            end
            plot!(p1.layout.grid[c_m,c_o],ylims=_ylims)
        end
    end
    savefig(p1,pth*".svg")
    return nothing
end

function permute_order1_coeffs(coeffs,model,o)
    return [
        coeffs[1:N₋(model)];
        permutedims(reshape(coeffs[N₋(model)+1:end-N₊(model)],n_phases(model),o,:),[2,1,3])[:];
        coeffs[end-N₊(model)+1:end];
    ]
end
# function make_plots!(orders,models,truth,truth_ci,err_data,errs_ci,pth,_xlims,_ylims,_xticks,_yticks,_phases)

#     k = keys(models[minimum(keys(models))])

#     p1 = plot(layout = (length(k),length(orders)))
#     # q1 = plot(layout = (length(k),length(orders)))

#     # what a truely terrible piece of code!
#     for (c_o,o) in enumerate(orders)
#         for (c_m,m) in enumerate(sort([i for i in k]))
#             for i in _phases
#                 # @eval pm = sum($(Symbol(m,"_π",o))[[1:2;end-1:end]])
#                 # std_plot!(p1.layout.grid[c_m,c_o],x->models[o][m]["pdf"](x,i),color=i)
#                 x_vals = range(_xlims[1],_xlims[2],length=200)
#                 y_vals = truth.pdf.(x_vals,i)
#                 ci = truth_ci.pdf.(x_vals,i)
#                 lwr = y_vals-[ci[n][1] for n in 1:length(ci)]
#                 upr = [ci[n][2] for n in 1:length(ci)]-y_vals
#                 plot!(p1.layout.grid[c_m,c_o],x_vals,y_vals,ribbon=(lwr,upr),color=i+2,label=false,alpha=0.5)
#                 std_plot!(p1.layout.grid[c_m,c_o],x->truth.pdf(x,i),color=i+2)
#                 # std_plot!(p1.layout.grid[c_m,c_o],x->truth_ci.pdf(x,i)[2],linestyle=:dot,color=i+2)

#                 std_plot!(q1.layout.grid[c_m,c_o],x->models[o][m]["cdf"](x,i),color=i)
#                 std_plot!(q1.layout.grid[c_m,c_o],x->truth.cdf(x,i),color=i+2)
#                 y_vals = truth.cdf.(x_vals,i)
#                 ci = truth_ci.cdf.(x_vals,i)
#                 lwr = y_vals-[ci[n][1] for n in 1:length(ci)]
#                 upr = [ci[n][2] for n in 1:length(ci)]-y_vals
#                 plot!(q1.layout.grid[c_m,c_o],x_vals,y_vals,ribbon=(lwr,upr),color=i+2,label=false,alpha=0.5)
#                 # std_plot!(q1.layout.grid[c_m,c_o],x->truth_ci.cdf(x,i)[1],linestyle=:dot,color=i+2)
#                 # std_plot!(q1.layout.grid[c_m,c_o],x->truth_ci.cdf(x,i)[2],linestyle=:dot,color=i+2)
#             end

#             if c_o==1
#                 plot!(p1.layout.grid[c_m,c_o],ylabel=uppercase(string(m)[1:end-1]))
#                 # plot!(q1.layout.grid[c_m,c_o],ylabel=uppercase(string(m)[1:end-1]))

#                 plot!(p1.layout.grid[c_m,c_o],yticks=_yticks[1])
#                 # plot!(q1.layout.grid[c_m,c_o],yticks=_yticks[2])
#             end
#             if c_m==1
#                 plot!(p1.layout.grid[c_m,c_o],title=string(o))
#                 # plot!(q1.layout.grid[c_m,c_o],title=string(o))
#             end
#             if c_m==4
#                 plot!(p1.layout.grid[c_m,c_o],xticks=_xticks[1])
#                 # plot!(q1.layout.grid[c_m,c_o],xticks=_xticks[2])
#             end
#             plot!(p1.layout.grid[c_m,c_o],ylims=_ylims[1])
#             # plot!(q1.layout.grid[c_m,c_o],ylims=_ylims[2])
#         end
#     end
#     file = pth*"/figs"
#     savefig(p1,file*"/pdfs.svg")
#     savefig(q1,file*"/cdfs.svg")
#     titles = [
#         "L¹ error between the true CDF and approximations";
#         "KS error between the true CDF and approximations";
#         "L¹ error between the true PDF and approximations";
#         "L² error between the true PDF and approximations";
#     ]
#     err_names = (:L1_cdf_errors,:ks_errors,:L1_pdf_errors,:L2_pdf_errors)
#     for (c_err,err) in enumerate(err_data)
#         q = plot()
#         for (c_m,m) in enumerate(names(err))
#             idx = [!(ismissing(err[i,m])||ismissing(errs_ci[1][c_err][i,m])||ismissing(errs_ci[2][c_err][i,m])) for i in 1:length(orders)]
#             lwr = err[idx,m]-errs_ci[1][c_err][idx,m]
#             upr = errs_ci[2][c_err][idx,m]-err[idx,m]
#             plot!(q,orders[idx],err[idx,m],ribbon=(lwr,upr),label=m,color=c_m)
#             # idx = [!ismissing(errs_ci[1][c_err][i,m]) for i in 1:length(orders)]
#             # plot!(q,orders[idx],errs_ci[1][c_err][idx,m],label=false,linestyle=:dot,color=c_m)
#             # idx = [!ismissing(errs_ci[2][c_err][i,m]) for i in 1:length(orders)]
#             # plot!(q,orders[idx],errs_ci[2][c_err][idx,m],label=false,linestyle=:dot,color=c_m)
#         end
#         plot!(q,xlabel="Order", ylabel="log₁₀ Error",
#             title=titles[c_err],
#             legend=:outertopright)
#         #display(q)
#         # #display((@__DIR__)*)
#         file = pth*"/figs/"*string(err_names[c_err])
#         savefig(q,file*".svg")
#     end
#     return p1, q1
# end
function ci_from_samples(samples,p)
    ci_lwr = Array{Float64,2}(undef,size(samples)[2:end]...)
    ci_upr = Array{Float64,2}(undef,size(samples)[2:end]...)
    for g in 1:size(samples,2)
        for i in 1:size(samples,3)
            ci_lwr[g,i], ci_upr[g,i] = quantile!(samples[:,g,i],p)
        end
    end
    return ci_lwr, ci_upr
end
function cdf_bootstrap_ci(sims,n_boot,grid,_phases,p,rng)
    n_sims = length(sims.φ)
    cdfs = Array{Float64,3}(undef,n_boot,length(grid),length(_phases))
    for n in 1:n_boot
        idx = rand(rng,1:n_sims,n_sims)
        boot_sample = Simulation(sims.t[idx],sims.φ[idx],sims.X[idx],sims.n[idx],sims.model)
        cdfs[n,:,:] = cdf(boot_sample).(grid[:],_phases[:]')
    end
    ci_lwr, ci_upr = ci_from_samples(cdfs,p)
    return ci_lwr, ci_upr
end
function cell_probs_bootstrap_ci(sims,n_boot,nodes,_phases,p,rng)
    n_sims = length(sims.φ)
    centres = (nodes[1:end-1]+nodes[2:end])/2.0
    cell_probs_evald = Array{Float64,3}(undef,n_boot,length(centres),length(_phases))
    for n in 1:n_boot
        idx = rand(rng,1:n_sims,n_sims)
        boot_sample = Simulation(sims.t[idx],sims.φ[idx],sims.X[idx],sims.n[idx],sims.model)
        cell_probs_evald[n,:,:] = cell_probs(boot_sample,nodes).(centres[:],_phases[:]')
    end
    ci_lwr, ci_upr = ci_from_samples(cell_probs_evald,p)
    return ci_lwr, ci_upr
end
function point_mass_boostrap_ci(sims,n_boot,p,rng)
    n_sims = length(sims.φ)
    pm_samps = Array{Float64,3}(undef,n_boot,1,N₋(sims.model)+N₊(sims.model))
    for n in 1:n_boot
        idx = rand(rng,1:n_sims,n_sims)
        boot_sample = Simulation(sims.t[idx],sims.φ[idx],sims.X[idx],sims.n[idx],sims.model)
        pm_samps[n,1,:] = sim_point_masses(boot_sample)
    end
    ci_lwr, ci_upr = ci_from_samples(pm_samps,p)
    return ci_lwr, ci_upr
end

function bootstrap_cis(sims,n_boot,grid,nodes,_phases,p,rng)
    n_sims = length(sims.φ)
    cdfs = Array{Float64,3}(undef,n_boot,length(grid),length(_phases))
    centres = (nodes[1:end-1]+nodes[2:end])/2.0
    cell_probs_evald = Array{Float64,3}(undef,n_boot,length(centres),length(_phases))
    pm_samps = Array{Float64,3}(undef,n_boot,1,N₋(sims.model)+N₊(sims.model))
    for n in 1:n_boot
        idx = rand(rng,1:n_sims,n_sims)
        boot_sample = Simulation(sims.t[idx],sims.φ[idx],sims.X[idx],sims.n[idx],sims.model)
        cdfs[n,:,:] = cdf(boot_sample).(grid[:],_phases[:]')
        cell_probs_evald[n,:,:] = cell_probs(boot_sample,nodes).(centres[:],_phases[:]')
        pm_samps[n,1,:] = sim_point_masses(boot_sample)
    end
    cdf_ci_lwr, cdf_ci_upr = ci_from_samples(cdfs,p)
    cell_probs_ci_lwr, cell_probs_ci_upr = ci_from_samples(cell_probs_evald,p)
    pm_ci_lwr, pm_ci_upr = ci_from_samples(pm_samps,p)
    return Dict{String,Any}(
        "bootstrap_samples" => Dict("cdf"=>cdfs,"cell_probs"=>cell_probs_evald,"boundary_masses"=>pm_samps),
        "cdf_ci_lwr" => cdf_ci_lwr,
        "cdf_ci_upr" => cdf_ci_upr,
        "cell_probs_ci_lwr" => cell_probs_ci_lwr,
        "cell_probs_ci_upr" => cell_probs_ci_upr,
        "boundary_masses_ci_lwr" => pm_ci_lwr[:],
        "boundary_masses_ci_upr" => pm_ci_upr[:],
    )
end

function write_models(models,pth)
    orders = sort([l for l in keys(models)])
    k = sort([k for k in keys(first(models)[2])])
    for (c_o,o) in enumerate(orders)
        for (c_m,m) in enumerate(k)
            temp_path = mkpath(pth*"/data/models/order_"*string(o)*"/model_"*string(m)*"/")
            for var in keys(models[o][m])
                try
                    JSON.write(temp_path*"/"*var*".json",JSON.json(models[o][m][var]))
                catch
                    JSON.write(temp_path*"/"*var*".json",JSON.json("failed to write "*var))
                end
            end
        end
    end
    return nothing
end
function write_truth(truth,pth)
    temp_path = mkpath(pth*"/data/sims/")
    for var in keys(truth)
        try
            JSON.write(temp_path*"/"*var*".json",JSON.json(truth[var]))
        catch
            JSON.write(temp_path*"/"*var*".json",JSON.json("failed to write "*var))
        end
    end
    return nothing
end

function read_directory_of_JSONs_to_dict(pth)

end
# function bootstrap_errs_ci(sims,orders,models,n_boot,n_evals,p,rng)
#     n_sims = length(sims.φ)
#     boot_errors = Array{Any,1}(undef,n_boot)
#     for n in 1:n_boot
#         idx = rand(rng,1:n_sims,n_sims)
#         boot_sample = Simulation(sims.t[idx],sims.φ[idx],sims.X[idx],sims.n[idx],sims.model)
#         truth =  (pdf=(x,i)->missing, cdf=cdf(boot_sample))
#         boot_errors[n] = compute_errors(orders,models,truth,n_evals)
#     end
#     k = keys(models[minimum(keys(models))])

#     L1_cdf_errors_lwr = DataFrame(DG = Vector{Union{Missing, Float64}}(undef,0), DG_Limit = Vector{Union{Missing, Float64}}(undef,0), Order_1 = Vector{Union{Missing, Float64}}(undef,0), QBDRAP = Vector{Union{Missing, Float64}}(undef,0))
#     ks_errors_lwr = DataFrame(DG = Vector{Union{Missing, Float64}}(undef,0), DG_Limit = Vector{Union{Missing, Float64}}(undef,0), Order_1 = Vector{Union{Missing, Float64}}(undef,0), QBDRAP = Vector{Union{Missing, Float64}}(undef,0))
#     L1_pdf_errors_lwr = DataFrame(DG = Vector{Union{Missing, Float64}}(undef,0), DG_Limit = Vector{Union{Missing, Float64}}(undef,0), Order_1 = Vector{Union{Missing, Float64}}(undef,0), QBDRAP = Vector{Union{Missing, Float64}}(undef,0))
#     L2_pdf_errors_lwr = DataFrame(DG = Vector{Union{Missing, Float64}}(undef,0), DG_Limit = Vector{Union{Missing, Float64}}(undef,0), Order_1 = Vector{Union{Missing, Float64}}(undef,0), QBDRAP = Vector{Union{Missing, Float64}}(undef,0))
#     L1_cdf_errors_upr = DataFrame(DG = Vector{Union{Missing, Float64}}(undef,0), DG_Limit = Vector{Union{Missing, Float64}}(undef,0), Order_1 = Vector{Union{Missing, Float64}}(undef,0), QBDRAP = Vector{Union{Missing, Float64}}(undef,0))
#     ks_errors_upr = DataFrame(DG = Vector{Union{Missing, Float64}}(undef,0), DG_Limit = Vector{Union{Missing, Float64}}(undef,0), Order_1 = Vector{Union{Missing, Float64}}(undef,0), QBDRAP = Vector{Union{Missing, Float64}}(undef,0))
#     L1_pdf_errors_upr = DataFrame(DG = Vector{Union{Missing, Float64}}(undef,0), DG_Limit = Vector{Union{Missing, Float64}}(undef,0), Order_1 = Vector{Union{Missing, Float64}}(undef,0), QBDRAP = Vector{Union{Missing, Float64}}(undef,0))
#     L2_pdf_errors_upr = DataFrame(DG = Vector{Union{Missing, Float64}}(undef,0), DG_Limit = Vector{Union{Missing, Float64}}(undef,0), Order_1 = Vector{Union{Missing, Float64}}(undef,0), QBDRAP = Vector{Union{Missing, Float64}}(undef,0))
#     for (c_o,o) in enumerate(orders)
#         L1_cdf_error_row_lwr = Vector{Union{Missing, Float64}}(undef,length(k))
#         ks_error_row_lwr = Vector{Union{Missing, Float64}}(undef,length(k))
#         L1_pdf_error_row_lwr = Vector{Union{Missing, Float64}}(undef,length(k))
#         L2_pdf_error_row_lwr = Vector{Union{Missing, Float64}}(undef,length(k))
#         L1_cdf_error_row_upr = Vector{Union{Missing, Float64}}(undef,length(k))
#         ks_error_row_upr = Vector{Union{Missing, Float64}}(undef,length(k))
#         L1_pdf_error_row_upr = Vector{Union{Missing, Float64}}(undef,length(k))
#         L2_pdf_error_row_upr = Vector{Union{Missing, Float64}}(undef,length(k))
#         for (c_m,m) in enumerate(sort([i for i in k]))
#             L1_cdf_error_row_lwr[c_m] = quantile(skipmissing([boot_errors[n][1][c_o,c_m] for n in 1:n_boot]),p[1])
#             ks_error_row_lwr[c_m] = quantile(skipmissing([boot_errors[n][2][c_o,c_m] for n in 1:n_boot]),p[1])
#             L1_pdf_error_row_lwr[c_m] = missing#quantile(skipmissing([boot_errors[n][3][c_o,c_m] for n in 1:n_boot]),p[1])
#             L2_pdf_error_row_lwr[c_m] = missing#quantile(skipmissing([boot_errors[n][4][c_o,c_m] for n in 1:n_boot]),p[1])
#             L1_cdf_error_row_upr[c_m] = quantile(skipmissing([boot_errors[n][1][c_o,c_m] for n in 1:n_boot]),p[2])
#             ks_error_row_upr[c_m] = quantile(skipmissing([boot_errors[n][2][c_o,c_m] for n in 1:n_boot]),p[2])
#             L1_pdf_error_row_upr[c_m] = missing#quantile(skipmissing([boot_errors[n][3][c_o,c_m] for n in 1:n_boot]),p[2])
#             L2_pdf_error_row_upr[c_m] = missing#quantile(skipmissing([boot_errors[n][4][c_o,c_m] for n in 1:n_boot]),p[2])
#         end
#         push!(L1_cdf_errors_lwr,L1_cdf_error_row_lwr)
#         push!(ks_errors_lwr,ks_error_row_lwr)
#         push!(L1_pdf_errors_lwr,L1_pdf_error_row_lwr)
#         push!(L2_pdf_errors_lwr,L2_pdf_error_row_lwr)
#         push!(L1_cdf_errors_upr,L1_cdf_error_row_upr)
#         push!(ks_errors_upr,ks_error_row_upr)
#         push!(L1_pdf_errors_upr,L1_pdf_error_row_upr)
#         push!(L2_pdf_errors_upr,L2_pdf_error_row_upr)
#     end

#     return (L1_cdf_errors_lwr,ks_errors_lwr,L1_pdf_errors_lwr,L2_pdf_errors_lwr), (L1_cdf_errors_upr,ks_errors_upr,L1_pdf_errors_upr,L2_pdf_errors_upr)
# end


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

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

# ╔═╡ 6b473b2d-4326-46b4-af38-07b61de287fc
begin
	import Pkg
	Pkg.activate(mktempdir())
	Pkg.add([
		Pkg.PackageSpec(name="ImageIO", version="0.5"),
		Pkg.PackageSpec(name="ImageShow", version="0.2"),
		Pkg.PackageSpec(name="FileIO", version="1.6"),
		Pkg.PackageSpec(name="PNGFiles", version="0.3.6"),
		Pkg.PackageSpec(name="Colors", version="0.12"),
		Pkg.PackageSpec(name="ColorVectorSpace", version="0.8"),

		Pkg.PackageSpec(name="PlutoUI", version="0.7"),
		Pkg.PackageSpec(name="HypertextLiteral", version="0.5"),
		Pkg.PackageSpec(name="ForwardDiff", version="0.10")
	])

	using Colors, ColorVectorSpace, ImageShow, FileIO
	using PlutoUI
	using HypertextLiteral
	using LinearAlgebra
	using ForwardDiff
end

# ╔═╡ d49682ff-d529-4283-871b-f8ee50a4e6ee
filter!(LOAD_PATH) do path
	path != "@v#.#"
end;

# ╔═╡ 2e8c4a48-d535-44ac-a1f1-4cb26c4aece6
filter!(LOAD_PATH) do path
	path != "@v#.#"
end;

# ╔═╡ c09c8ba0-887e-11eb-07e3-71377ec0e708
html"""

Section 1.3

 Transformations & Autodiff

Lecture Video

"""

# ╔═╡ b7895bd2-7634-11eb-211e-ef876d23bd88
PlutoUI.TableOfContents(aside=true)

# ╔═╡ e6a09409-f262-453b-a434-bfd935306719
md"""
#### Intializing packages

_When running this notebook for the first time, this could take up to 15 minutes. Hang in there!_
"""

# ╔═╡ 58a520ca-763b-11eb-21f4-3f27aafbc498
md"""
Last time, recall we defined linear combinations of images.  Remember we
* scaled an image by  multiplying by a constant
* combined images by adding the colors in each pixel possibly saturating
In general if we perform both of these operations, we get a linear combination.
"""

# ╔═╡ 2cca0638-7635-11eb-3b60-db3fabe6f536
md"""
# 4.1 Functions in Math and Julia
"""

# ╔═╡ c8a3b5b4-76ac-11eb-14f0-abb7a33b104d
md"""
### 4.1.1 Univariate Functions
"""

# ╔═╡ db56bcda-76aa-11eb-2447-5d9076789244
md"""
In high school you learned about univariate functions e.g.
* $f₁(x)=x^2$
* $f₂(x)=\sin(x)$
* $f₃(x)=x^\alpha$

In Julia, functions can be written in short form, anonymous form, or long form.
"""

# ╔═╡ 539aeec8-76ab-11eb-32a3-95c6672a0ea9
# short form
f₁(x) = x^2 # subscript unicode:   \_1 + <tab>

# ╔═╡ 81a00b78-76ab-11eb-072a-6b96847c2ce4
f₁(5)

# ╔═╡ 2369fb18-76ab-11eb-1189-85309c8f925b
# anonymous form
x->sin(x)

# ╔═╡ 98498f84-76ab-11eb-23cf-857c776a9163
(x->sin(x))(π/2) # <=> g(π/2) where g = x->sin(x)

# ╔═╡ c6c860a6-76ab-11eb-1dec-1b2f179a0fa9
f₃(x, α=3) = x^α

# ╔═╡ f07fbc6c-76ab-11eb-3382-87c7d65b4078
 f₃(5)

# ╔═╡ f4fa8c1a-76ab-11eb-302d-bd410432e3cf
f₃(5, 2)

# ╔═╡ b3faf4d8-76ac-11eb-0be9-7dda3d37aba0
md"""
Keywords
"""

# ╔═╡ 71c074f0-76ac-11eb-2164-c36381212bff
f₄(x;α) = x^α

# ╔═╡ 87b99c8a-76ac-11eb-1c94-8f1ffe3be593
f₄(2; α=5)

# ╔═╡ 504076fc-76ac-11eb-30c3-bfa75c991cb2
md"""
See [Julia's function documentation](https://docs.julialang.org/en/v1/manual/functions/) for more.
"""

# ╔═╡ f1dd24d8-76ac-11eb-1de7-a763a1b95668
md"""
### 4.1.2 Automatic Differentiation of Univariates
"""

# ╔═╡ fe01da74-76ac-11eb-12e3-8320340b6139
md"""
Automatic differentiation is a key enabling technology for machine learning and so much of scientific computing.  It derives the answer in a manner that is interestingly different from the symbolic differentiation of elementary calculus classes and the numerical differentiation of simple finite differences.  See the video at the end of this lecture.
"""

# ╔═╡ d42aec08-76ad-11eb-361a-a1f2c90fd4ec
ForwardDiff.derivative(f₁, 5)

# ╔═╡ 85d96989-b9eb-454b-9ca7-6454f19de12e
ForwardDiff.derivative(x->f₃(x, 3), 5) == 3 * 5^2

# ╔═╡ 28cd454c-76ae-11eb-0d1e-a56995100d59
md"""
Notice the use of anonymous functions to fix the parameter α=3
"""

# ╔═╡ 38b51946-76ae-11eb-2c8a-e19b30bf42cb
md"""
In case you have forgotten what a derivative is, we remind you with a simple finite difference approximation:
"""

# ╔═╡ 632a1f8c-76ae-11eb-2088-15c3e3c0a210
begin
	md"""
	$(@bind e Slider(-6:-1, default=-1, show_value=true))
	"""
end

# ╔═╡ 8a99f186-76af-11eb-031b-f1c288993c7f
ϵ = 10.0^e

# ╔═╡ ca1dfb8a-76b0-11eb-1323-d100bdeedc2d
(sin(1+ϵ)-sin(1))/ϵ , cos(1), ForwardDiff.derivative(sin,1)

# ╔═╡ f7df6cda-76b1-11eb-11e4-8d0af0349651
md"""
### 4.1.3 Scalar Valued Multivariate Functions
"""

# ╔═╡ 63449b54-76b4-11eb-202f-3bda2f4cff4d
md"""
Sometimes we are interested in scalar valued functions of more than 1 variable.
This can be written in Julia as a function of many variables or a function of
a vector.
e.g. $f_5(x) = 5\sin(x_1*x_2) + 2x_2/4x_3$
"""

# ╔═╡ 8c6b0236-76b4-11eb-2acf-91da23bedf0e
begin
	f₅(v) = 5sin(v[1]*v[2]) + 2*v[2]/4v[3]
	f₅(x, y, z) = 5sin(x*y) + 2*y/4z
end

# ╔═╡ a397d526-76b5-11eb-3cce-4374e33324d1
f₅(1, 2, 3), f₅([1, 2, 3])

# ╔═╡ 4a57d898-76b6-11eb-15ea-7be43393922c
md"""
Better yet if you must write it the two ways ( you probably won't need to, but if you must), **don't copy code**, reuse code so if it changes in one place it changes everywhere.
"""

# ╔═╡ bf23ab30-76b5-11eb-1adb-3d74a52cddfd
begin
	f₆(x, y, z) = 5sin(x*y) + 2*y/4z
	f₆(v) = f₆(v[1],v[2],v[3])
end

# ╔═╡ d5d4ac48-76b6-11eb-1687-ed853c2db7c9
f₆(1,2,3), f₆([1,2,3])

# ╔═╡ 89b2d570-76ba-11eb-0389-813bbb33efea
md"""
There's one other julia idea that is a trick to make vector code more readable. If you give a tuple argument, the function works directly on vectors but is defined with readable letters.
"""

# ╔═╡ a8c28578-76ba-11eb-3f3f-af35ff0b6c74
f₇((x,y,z)) = 5sin(x*y) + 2*y/4z # more readable then 5sin(v[1]*v[2]) + 2*v[2]/4v[3]

# ╔═╡ d9e07084-76ba-11eb-18ac-c58b1bc972ba
f₇([1,2,3]) # this works with vector arguments, but not scalars (f₇(1,2,3) errros)

# ╔═╡ 42172fb6-76b7-11eb-0a11-b7e10c6881f5
md"""
You can see that the functions $f_5$  and $f_6$ has two julia methods, one with one variable, and one with three variables.
"""

# ╔═╡ 57b0cd3c-76b7-11eb-0ece-810f3a8ede00
methods(f₅)

# ╔═╡ 6d411cea-76b9-11eb-061b-87d472bc3bdd
md"""
### 4.1.4 Automatic Differentiation: Scalar valued multivariate functions
"""

# ╔═╡ bc2c6afc-76b7-11eb-0631-51f83cd73205
md"""
In many applications, including machine learning, one needs to take derivatives of the function in every argument direction.  This is known as the *gradient*.  Automatic differentiation works again:
"""

# ╔═╡ ef06cfd8-76b7-11eb-1530-1fcd7e5c5992
ForwardDiff.gradient(f₅, [1,2,3])

# ╔═╡ 051db7a0-76b8-11eb-14c7-531f42ef60b8
md"""
Remember
$f_5(x) = 5\sin(x_1*x_2) + 2x_2/4x_3$
"""

# ╔═╡ 5f1afd24-76b8-11eb-36ab-9bbb3d73b930
md"""
One can check numerically by adding a small change to each of the arguments.m
"""

# ╔═╡ 2705bf34-76b8-11eb-3aaa-d363085784ff
begin
	∂f₅∂x = (f₅(1 + ϵ, 2, 3) - f₅(1, 2, 3)) / ϵ
	∂f₅∂y = (f₅(1, 2 + ϵ, 3) - f₅(1, 2, 3)) / ϵ
	∂f₅∂z = (f₅(1, 2, 3 + ϵ) - f₅(1, 2, 3)) / ϵ
	∇f = [∂f₅∂x, ∂f₅∂y, ∂f₅∂z]
end

# ╔═╡ dfb9d74c-76b8-11eb-24ff-e521f1294a6f
md"""
Whether you are an expert at multivariable calculus, or you have never seen this before, I hope seeing it numerically makes the idea intuitive.
"""

# ╔═╡ 1049f458-76b9-11eb-1d2d-af0b22480121
md"""
**Important Remark**: In machine learning, and other optimization contexts, we want to minimize a scalar function of many parameters known as a "loss function."  Following the negative gradient is a standard technique for minimizing functions especially when there are many variables.  When there are only a few variables, there are better techniques.
"""

# ╔═╡ a0afe3ae-76b9-11eb-2301-cde7260ddd7f
md"""
### 4.1.5. Transformations: Vector Valued Multivariate Functions
"""

# ╔═╡ ac1ab224-76bb-11eb-13cb-0bd44bea1042
md"""
While scalar functions might technically be called a transformation, it is more common to use the term when both the input and output are multidimensional.
"""

# ╔═╡ bcf92688-76b9-11eb-30fb-1f320a65f45a
md"""
Let us consider a few functions that take in a vector of size 2 and returns a vector of size 2.
"""

# ╔═╡ d364f91a-76b9-11eb-1807-75e733940d53
begin
	idy((x,y)) = [x,y]
	lin1((x,y)) =  [2x + 3y, -5x+4x]
	scalex(α) = ((x, y), ) -> (α * x, y)
	scaley(α) = ((x, y), ) -> (x, α * y)
	rot(θ) = ((x, y), ) -> [cos(θ)*x + sin(θ)*y, -sin(θ)*x + cos(θ)*y]
	shear(α) = ((x, y), ) -> [x + α * y, y]
	genlin(a, b, c, d) = ((x, y), ) -> [a*x + b*y; c*x + d*y]
end

# ╔═╡ 98936c7f-a628-4338-8b53-529a09199ec2
md"""
	Testing...
"""

# ╔═╡ 28027825-9dd7-40ff-a928-2687298b445e
idy([1, 2])

# ╔═╡ c85584e2-0a41-4923-a886-a7c43250565b
lin1([1, 2])

# ╔═╡ 392320b1-f8c6-4831-9a74-047782d6f2f4
scalex(π)([1, 2]) == (π * 1, 2)

# ╔═╡ c3b1012d-0ee2-450e-8f56-f5998a848ede
scaley(π)([1, 2]) == (1, π * 2)

# ╔═╡ f25c6308-76b9-11eb-3563-1f0ef4cdf86a
rot(π/2)([4, 5]) == [cos(π/2)*4 + sin(π/2)*5, -sin(π/2)*4 + cos(π/2)*5]

# ╔═╡ 0fbba269-a289-4cea-8fe9-79277bee8190
shear(π/2)([4, 5]) == [4 + π/2 * 5, 5]

# ╔═╡ 127e28ff-c25b-400c-9151-6d9ad03f9602
genlin(1, 2, 3, 4)([4, 5]) == [1*4 + 2*5; 3*4 + 4*5]

# ╔═╡ c9a148f0-76bb-11eb-0778-9d3e84369a19
md"""
We bet you have noticed that these functions could all have been defined with matrices. Indeed the general case can be written
"""

# ╔═╡ db4bc328-76bb-11eb-28dc-eb9df8892d01
# [a b;c d] * [x,y]

# ╔═╡ 89f0bc54-76bb-11eb-271b-3190b4d8cbc0
md"""
or in math `` \begin{pmatrix} a & b \\ c & d \end{pmatrix}
\begin{pmatrix} x \\ y \end{pmatrix}`` .
"""

# ╔═╡ f70f7ea8-76b9-11eb-3bd7-87d40a2861b1
md"""
By contrast here are a few fun functions that can not be written as matrix times
vector.  What characterizes the matrix ones from the non-matrix ones?
"""

# ╔═╡ 78176284-76bc-11eb-3045-f584127f58b9
begin
	function warp(α)
		((x,y),)  -> begin
			r = √(x^2+y^2)
			θ=α*r
			rot(θ)([x,y])
		end
	end

	rθ(x) = (norm(x), atan(x[2],x[1])) # maybe vectors are more readable here?

	xy((r,θ)) = (r*cos(θ), r*sin(θ))
end

# ╔═╡ bf28c388-76bd-11eb-08a7-af2671218017
md"""
This may be a little fancy, but we see that warp is a rotation, but the
rotation depends on the point where it is applied.
"""

# ╔═╡ 5655d2a6-76bd-11eb-3042-5b2dd3f6f44e
begin
	warp₂(α,x,y) = rot(α*√(x^2+y^2))
	warp₂(α) = ((x,y),) -> warp₂(α,x,y)([x,y])
end

# ╔═╡ 852592d6-76bd-11eb-1265-5f200e39113d
warp(1)([5,6])

# ╔═╡ 8e36f4a2-76bd-11eb-2fda-9d1424752812
warp₂(1.0)([5.0,6.0])

# ╔═╡ 09ed6d38-76be-11eb-255b-3fbf76c21097
md"""
### 4.1.6 Automatic Differentiation of Transformations
"""

# ╔═╡ 9786e2be-76be-11eb-3755-b5669c37aa64
ForwardDiff.jacobian( warp(3.0), [4,5] )

# ╔═╡ 963694d6-76be-11eb-1b27-d5d063964d24
md"""
What is this thing?
"""

# ╔═╡ b78ef2fe-76be-11eb-1f55-3d0874b298e8
begin
	∂w∂x = (warp(3.0)([4+ϵ, 5]) -   warp(3.0)([4,5]))/ϵ # This is a vector, right?
	∂w∂y = (warp(3.0)([4,   5+ϵ]) - warp(3.0)([4,5]))/ϵ # This too
	[∂w∂x ∂w∂y]
end

# ╔═╡ ad728ee6-7639-11eb-0b23-c37f1366fb4e
md"""
## 4.2 But what is a transformation, really?
You have very likely learned how to multiply matrices times vectors.  I'll bet you think of a matrix as a table of numbers, and a vector as a column of numbers, and if you are well practiced, you know just when to multiply and just when to add.
Congratulations, you now can do what computers excel at.

"""

# ╔═╡ 4d4e6b32-763b-11eb-3021-8bc61ac07eea
md"""
Matrices are often thought of as containers of numbers in a rectangular array, and hence one thinks of manipulating these tables like a spreadsheet, but actually the deeper meaning is that it is a transformation.
"""

# ╔═╡ ce55beee-7643-11eb-04bc-b517703facff
md"""
α= $(@bind α Slider(.1:.1:3, show_value=true))
"""

# ╔═╡ 005ca75a-7622-11eb-2ba4-9f450e71df1f
let

range = -1.5:.1:1.5
md"""
This is a "scrubbable matrix" -- click on the number and drag to change.

``(``
 $(@bind a Scrubbable( range; default=1.0))
 $(@bind b Scrubbable( range; default=0.0))
``)``

``(``
$(@bind c Scrubbable(range; default=0.0 ))
$(@bind d Scrubbable(range; default=1.0))
``)``

	**Re-run this cell to reset to identity transformation**
"""
end

# ╔═╡ ed3caab2-76bf-11eb-2544-21e8181adef5
#T = shear(1) # Pick a transformation
T = genlin(a,b,c,d)

# ╔═╡ 2efaa336-7630-11eb-0c17-a7d4a0141dac
md"""
center zoom = $(@bind z Slider(.1:.1:3, show_value=true, default=1))
"""

# ╔═╡ 488d732c-7631-11eb-38c3-e9a9165d5606
md"""
top left zoom =	$(@bind f Slider(.1:1:3, show_value=true, default=1))
"""

# ╔═╡ 60532aa0-740c-11eb-0402-af8ff117f042
md"Show grid lines $(@bind show_grid CheckBox(default=true))"

# ╔═╡ f085296d-48b1-4db6-bb87-db863bb54049
A = [
	a b
	c d
	]

# ╔═╡ d1757b2c-7400-11eb-1406-d937294d5388
md"**_Det(A)_ = $a * $d - $c * $b =  $(det(A))**"

# ╔═╡ 5227afd0-7641-11eb-0065-918cb8538d55
md"""
We never seem to see this in linear algebra classes do we?

Check out
[Linear Map Wikipedia](https://en.wikipedia.org/wiki/Linear_map)

[Transformation Matrix Wikipedia](https://en.wikipedia.org/wiki/Transformation_matrix)
"""

# ╔═╡ 2835e33a-7642-11eb-33fd-79fb8ad27fa7
md"""
Geometry of determinant, how areas scale.
"""

# ╔═╡ 4c93d784-763d-11eb-1f48-81d4d45d5ce0
md"""
## Why are we doing this backwards?
"""

# ╔═╡ 559660cc-763d-11eb-1ed8-017b93ed4ecb
md"""
Computer Science
Solving 2 equations in 2 unknowns, and higher dimensional analogs.

The top 500 supercomputers, and how many equations in how many unknowns are being solved today.
"""

# ╔═╡ a66eb6fe-76b3-11eb-1d50-659ec2bf7c44
md"""
### 4.4 Automatic Differentiation in 10 mins(ok 11)
"""

# ╔═╡ b9dba026-76b3-11eb-1bfb-ffe9c43ced5d
html"""

"""

# ╔═╡ 62b28c02-763a-11eb-1418-c1e30555b1fa
md"""
What about more than 3 dimensions?

Oh I've heard it all,
* there is no such thing as more than 3 dimensions.
* The fourth dimension is time, isn't it?
* I have trouble visualizing 3 dimensions, let alone contemplating 4.

... but mathematicians have no trouble with n spatial dimensions, and they do not (usually) worry about silly questions like whether these exist.
"""

# ╔═╡ c536dafb-4206-4689-ad6d-6935385d8fdf
md"""
## Appendix
"""

# ╔═╡ fb509fb4-9608-421d-9c40-a4375f459b3f
det_A = det(A)

# ╔═╡ 40655bcc-6d1e-4d1e-9726-41eab98d8472
img_sources = [
	"https://user-images.githubusercontent.com/6933510/108605549-fb28e180-73b4-11eb-8520-7e29db0cc965.png" => "Corgis",
	"https://user-images.githubusercontent.com/6933510/108883855-39690f80-7606-11eb-8eb1-e595c6c8d829.png" => "Arrows",
	"https://images.squarespace-cdn.com/content/v1/5cb62a904d546e33119fa495/1589302981165-HHQ2A4JI07C43294HVPD/ke17ZwdGBToddI8pDm48kA7bHnZXCqgRu4g0_U7hbNpZw-zPPgdn4jUwVcJE1ZvWQUxwkmyExglNqGp0IvTJZamWLI2zvYWH8K3-s_4yszcp2ryTI0HqTOaaUohrI8PISCdr-3EAHMyS8K84wLA7X0UZoBreocI4zSJRMe1GOxcKMshLAGzx4R3EDFOm1kBS/fluffy+corgi?format=2500w" => "Long Corgi"
]

# ╔═╡ c0c90fec-0e55-4be3-8ea2-88b8705ee258
md"""
#### Choose an image:

$(@bind img_source Select(img_sources))
"""

# ╔═╡ 4fcb4ac1-1ad1-406e-8776-4675c0fdbb43
img_original = load(download(img_source));

# ╔═╡ 52a8009e-761c-11eb-2dc9-dbccdc5e7886
typeof(img_original)

# ╔═╡ b754bae2-762f-11eb-1c6a-01251495a9bb
begin
	white(c::RGB) = RGB(1,1,1)
	white(c::RGBA) = RGBA(1,1,1,0.75)
end

# ╔═╡ 7d0096ad-d89a-4ade-9679-6ee95f7d2044
function trygetpixel(img::AbstractMatrix, x::Float64, y::Float64)
	rows, cols = size(img)

	"The linear map [-1,1] ↦ [0,1]"
	f = t -> (t - -1.0)/(1.0 - -1.0)

	i = floor(Int, rows *  f(-y) / z)
	j = floor(Int, cols *  f(x * (rows / cols))  / z)

	if 1 < i ≤ rows && 1 < j ≤ cols
		img[i,j]
	else
		white(img[1,1])

	end
end

# ╔═╡ 83d45d42-7406-11eb-2a9c-e75efe62b12c
function with_gridlines(img::Array{<:Any,2}; n=16)

	sep_i = size(img, 1) ÷ n
	sep_j = size(img, 2) ÷ n

	result = copy(img)
	# stroke = zero(eltype(img))#RGBA(RGB(1,1,1), 0.75)

	stroke = RGBA(1, 1, 1, 0.75)

	result[1:sep_i:end, :] .= stroke
	result[:, 1:sep_j:end] .= stroke

	# a second time, to create a line 2 pixels wide
	result[2:sep_i:end, :] .= stroke
	result[:, 2:sep_j:end] .= stroke

	 result[  sep_i * (n ÷2) .+ [1,2]    , :] .= RGBA(0,1,0,1)
	result[ : ,  sep_j * (n ÷2) .+ [1,2]    ,] .= RGBA(1,0,0,1)
	return result
end

# ╔═╡ 55898e88-36a0-4f49-897f-e0850bd2b0df
img = if show_grid
	with_gridlines(img_original)
else
	img_original
end;

# ╔═╡ 8e0505be-359b-4459-9de3-f87ec7b60c23
[
	if det_A == 0
		RGB(1.0, 1.0, 1.0)
	else

		 # in_x, in_y = A \ [out_x, out_y]
         # in_x, in_y = xy( [out_x, out_y] )
		in_x, in_y =  T([out_x, out_y])
		trygetpixel(img, in_x, in_y)
	end

	for out_y in LinRange(f, -f, 500),
		out_x in LinRange(-f, f, 500)
]

# ╔═╡ Cell order:
# ╟─c09c8ba0-887e-11eb-07e3-71377ec0e708
# ╟─b7895bd2-7634-11eb-211e-ef876d23bd88
# ╟─e6a09409-f262-453b-a434-bfd935306719
# ╠═6b473b2d-4326-46b4-af38-07b61de287fc
# ╟─d49682ff-d529-4283-871b-f8ee50a4e6ee
# ╟─58a520ca-763b-11eb-21f4-3f27aafbc498
# ╟─2cca0638-7635-11eb-3b60-db3fabe6f536
# ╟─c8a3b5b4-76ac-11eb-14f0-abb7a33b104d
# ╟─db56bcda-76aa-11eb-2447-5d9076789244
# ╠═539aeec8-76ab-11eb-32a3-95c6672a0ea9
# ╠═81a00b78-76ab-11eb-072a-6b96847c2ce4
# ╠═2369fb18-76ab-11eb-1189-85309c8f925b
# ╠═98498f84-76ab-11eb-23cf-857c776a9163
# ╠═c6c860a6-76ab-11eb-1dec-1b2f179a0fa9
# ╠═f07fbc6c-76ab-11eb-3382-87c7d65b4078
# ╠═f4fa8c1a-76ab-11eb-302d-bd410432e3cf
# ╟─b3faf4d8-76ac-11eb-0be9-7dda3d37aba0
# ╠═71c074f0-76ac-11eb-2164-c36381212bff
# ╠═87b99c8a-76ac-11eb-1c94-8f1ffe3be593
# ╟─504076fc-76ac-11eb-30c3-bfa75c991cb2
# ╟─f1dd24d8-76ac-11eb-1de7-a763a1b95668
# ╟─fe01da74-76ac-11eb-12e3-8320340b6139
# ╠═d42aec08-76ad-11eb-361a-a1f2c90fd4ec
# ╠═85d96989-b9eb-454b-9ca7-6454f19de12e
# ╟─28cd454c-76ae-11eb-0d1e-a56995100d59
# ╟─38b51946-76ae-11eb-2c8a-e19b30bf42cb
# ╟─632a1f8c-76ae-11eb-2088-15c3e3c0a210
# ╟─8a99f186-76af-11eb-031b-f1c288993c7f
# ╠═ca1dfb8a-76b0-11eb-1323-d100bdeedc2d
# ╟─f7df6cda-76b1-11eb-11e4-8d0af0349651
# ╟─63449b54-76b4-11eb-202f-3bda2f4cff4d
# ╠═8c6b0236-76b4-11eb-2acf-91da23bedf0e
# ╠═a397d526-76b5-11eb-3cce-4374e33324d1
# ╟─4a57d898-76b6-11eb-15ea-7be43393922c
# ╠═bf23ab30-76b5-11eb-1adb-3d74a52cddfd
# ╠═d5d4ac48-76b6-11eb-1687-ed853c2db7c9
# ╟─89b2d570-76ba-11eb-0389-813bbb33efea
# ╠═a8c28578-76ba-11eb-3f3f-af35ff0b6c74
# ╠═d9e07084-76ba-11eb-18ac-c58b1bc972ba
# ╟─42172fb6-76b7-11eb-0a11-b7e10c6881f5
# ╠═57b0cd3c-76b7-11eb-0ece-810f3a8ede00
# ╟─6d411cea-76b9-11eb-061b-87d472bc3bdd
# ╟─bc2c6afc-76b7-11eb-0631-51f83cd73205
# ╠═ef06cfd8-76b7-11eb-1530-1fcd7e5c5992
# ╟─051db7a0-76b8-11eb-14c7-531f42ef60b8
# ╟─5f1afd24-76b8-11eb-36ab-9bbb3d73b930
# ╠═2705bf34-76b8-11eb-3aaa-d363085784ff
# ╟─dfb9d74c-76b8-11eb-24ff-e521f1294a6f
# ╟─1049f458-76b9-11eb-1d2d-af0b22480121
# ╟─a0afe3ae-76b9-11eb-2301-cde7260ddd7f
# ╟─ac1ab224-76bb-11eb-13cb-0bd44bea1042
# ╟─bcf92688-76b9-11eb-30fb-1f320a65f45a
# ╠═d364f91a-76b9-11eb-1807-75e733940d53
# ╟─98936c7f-a628-4338-8b53-529a09199ec2
# ╠═28027825-9dd7-40ff-a928-2687298b445e
# ╠═c85584e2-0a41-4923-a886-a7c43250565b
# ╠═392320b1-f8c6-4831-9a74-047782d6f2f4
# ╠═c3b1012d-0ee2-450e-8f56-f5998a848ede
# ╠═f25c6308-76b9-11eb-3563-1f0ef4cdf86a
# ╠═0fbba269-a289-4cea-8fe9-79277bee8190
# ╠═127e28ff-c25b-400c-9151-6d9ad03f9602
# ╟─c9a148f0-76bb-11eb-0778-9d3e84369a19
# ╠═db4bc328-76bb-11eb-28dc-eb9df8892d01
# ╟─89f0bc54-76bb-11eb-271b-3190b4d8cbc0
# ╟─f70f7ea8-76b9-11eb-3bd7-87d40a2861b1
# ╠═78176284-76bc-11eb-3045-f584127f58b9
# ╟─bf28c388-76bd-11eb-08a7-af2671218017
# ╠═5655d2a6-76bd-11eb-3042-5b2dd3f6f44e
# ╠═852592d6-76bd-11eb-1265-5f200e39113d
# ╠═8e36f4a2-76bd-11eb-2fda-9d1424752812
# ╟─09ed6d38-76be-11eb-255b-3fbf76c21097
# ╠═9786e2be-76be-11eb-3755-b5669c37aa64
# ╟─963694d6-76be-11eb-1b27-d5d063964d24
# ╠═b78ef2fe-76be-11eb-1f55-3d0874b298e8
# ╟─ad728ee6-7639-11eb-0b23-c37f1366fb4e
# ╟─4d4e6b32-763b-11eb-3021-8bc61ac07eea
# ╟─2e8c4a48-d535-44ac-a1f1-4cb26c4aece6
# ╟─c0c90fec-0e55-4be3-8ea2-88b8705ee258
# ╟─ce55beee-7643-11eb-04bc-b517703facff
# ╟─005ca75a-7622-11eb-2ba4-9f450e71df1f
# ╠═ed3caab2-76bf-11eb-2544-21e8181adef5
# ╟─2efaa336-7630-11eb-0c17-a7d4a0141dac
# ╟─488d732c-7631-11eb-38c3-e9a9165d5606
# ╟─60532aa0-740c-11eb-0402-af8ff117f042
# ╠═8e0505be-359b-4459-9de3-f87ec7b60c23
# ╟─f085296d-48b1-4db6-bb87-db863bb54049
# ╟─d1757b2c-7400-11eb-1406-d937294d5388
# ╟─5227afd0-7641-11eb-0065-918cb8538d55
# ╟─2835e33a-7642-11eb-33fd-79fb8ad27fa7
# ╟─4c93d784-763d-11eb-1f48-81d4d45d5ce0
# ╟─559660cc-763d-11eb-1ed8-017b93ed4ecb
# ╟─a66eb6fe-76b3-11eb-1d50-659ec2bf7c44
# ╟─b9dba026-76b3-11eb-1bfb-ffe9c43ced5d
# ╟─62b28c02-763a-11eb-1418-c1e30555b1fa
# ╟─c536dafb-4206-4689-ad6d-6935385d8fdf
# ╟─fb509fb4-9608-421d-9c40-a4375f459b3f
# ╟─40655bcc-6d1e-4d1e-9726-41eab98d8472
# ╠═4fcb4ac1-1ad1-406e-8776-4675c0fdbb43
# ╠═52a8009e-761c-11eb-2dc9-dbccdc5e7886
# ╠═55898e88-36a0-4f49-897f-e0850bd2b0df
# ╠═7d0096ad-d89a-4ade-9679-6ee95f7d2044
# ╠═b754bae2-762f-11eb-1c6a-01251495a9bb
# ╠═83d45d42-7406-11eb-2a9c-e75efe62b12c

ext: jl


text:
function reverse(mg::MetaDiGraph)
    rg = reverse(mg.graph)
    rvprops = copy(mg.vprops)
    reprops = empty(mg.eprops)
    rgprops = mg.gprops
    rweightfunction = mg.weightfunction
    rdefaultweight = mg.defaultweight
    rindices = copy(mg.metaindex)

    for (u, v) in keys(mg.eprops)
        reprops[(v, u)] = mg.eprops[(u, v)]
    end

    rmg =
        MetaGraph(rg, rvprops, reprops, rgprops, rweightfunction, rdefaultweight, rindices)

    return rmg
end


text:
using Printf

function conditional_examples()
    println("--------------------------------------")
    println("Conditional examples")

    println()
    println("We can use all the usual && || ! >= <= != operators")
    age = 12
    if age >= 5 && age <= 6
        println("You're in Kindergarten")
    elseif age >= 7 && age <= 13
        println("You're in Middle School")
    elseif age >= 14 && age <= 18
        println("You're in High School")
    else
        println("Stay Home")
    end

    println("We can use the @printf macro for evaluting (note that all macros start with the @-sign)")
    @printf("true || false = %s\n", true || false ? "true" : "false")
    @printf("!true = %s\n", !true ? "true" : "false")
    age = 55
    @printf("Can I drive : %s\n", age > 16 ? "true" : "false")
end

max_stars_repo_name: James-P-D/JuliaDump

max_issues_repo_name: James-P-D/JuliaDump

max_forks_repo_name: James-P-D/JuliaDump


text:
using Artifacts, ArtifactUtils

# update to newest version of PGLIB to create new artifact
# Example at ArtifactUtils.jl github
add_artifact!(
    "Artifacts.toml",
    "PGLib_opf",
    "https://github.com/power-grid-lib/pglib-opf/archive/refs/tags/v21.07.tar.gz",
    force=true,
)


text:
# poissonreg.jl
#
# Code for poisson regression.

using Distributions;

# (A, b, x) = GeneratePoissonData(m::Int64 = 200, n::Int64 = 20)
#
# Generates an m-by-n matrix A of measurement vectors a_i = A[i, :], where each
# a_i is uniform on the unit sphere. Chooses x uniformly from the sphere of
# radius sqrt(n), and draws b_i as
#
#   b_i ~ Poisson(exp(a_i' * x))
function GeneratePoissonData(m::Int64 = 200, n::Int64 = 20;
                             generate_gaussian::Bool = false,
                             condition_number::Float64 = 1.0)
  A = randn(m, n);
  if (!generate_gaussian)
    A = A ./ sqrt.(repeat(sum(A.^2, dims=2), outer = (1, n)));
  end
  if (condition_number > 1.0)
    # Make A have condition number given
    condition_vec = Vector(LinRange(1.0, condition_number, n));
    condition_vec = condition_vec / norm(condition_vec);
    A = A .* repeat(condition_vec', outer = (m, 1));
  end
  x = randn(n);
  x = sqrt(n) * x / norm(x);
  b = zeros(Int64, m);
  for ii = 1:m
    R = Poisson(exp(dot(A[ii, :], x)));
    b[ii] = rand(R);
  end
  return (A, b, x);
end

# PoissonObj(A::Matrix{Float64}, b::Vector{Int64}, x::Vector{Float64})
#
# Computes the (average) log loss for the poisson regression problem,
# where we recall that the log-likelihood for b ~ Poisson(lambda) is
#
#  p(b | lambda) = exp(-lambda) * lambda^b / b!
#
# For a data vector a, the parameter lambda = exp(a' * x), so we have log loss
# (ignoring the factorial term)
#
# l(x; (a, b)) = -b * dot(a, x) + exp(dot(a, x)).
function PoissonObj(A::Matrix{Float64}, b::Vector{Int64}, x::Vector{Float64})
  mm = size(A, 1);
  objective = -sum(b .* (A * x)) + sum(exp.(A*x)); # no need to add b!
  return objective / mm;
end

# g = FullPoissonGradient(A, b, x)
#
# Sets g to be the gradient of the poisson objective and returns it.
function FullPoissonGradient(A::Matrix{Float64}, b::Vector{Int64},
                             x::Vector{Float64})
  mm = size(A, 1);
  grad = (-A' * b + A' * exp.(A * x)) / mm;
  return grad;
end

# x = ProximalUpdatePoisson(a::Vector{Float64}, b::Float64,
#                              x_init::Vector, stepsize::Float64; accuracy)
#
# Sets x to minimize
#
#  -b * a' * x + exp(x' * a) + norm(x - x_init)^2 / (2 * stepsize).
#
# Returns the minimizing x. The accuracy parameter governs the
# accuracy of the solution: we recommend a line search.
function ProximalUpdatePoisson(a::Vector{Float64}, b::Int64,
                               x_init::Vector{Float64}, stepsize::Float64;
                               eps_accuracy::Float64 = 1e-10)
  # Replace this code to return the correct update
  return x_init;
end

# x = TruncatedUpdatePoisson(a::Vector{Float64}, b::Int64,
#                            x_init::Vector{Float64}, stepsize::Float64)
#
# Let l(x; (a, b)) = -b * a' * x + exp(a' * x). Sets x to minimize the
# lower-bounded approximation to F at x_init, that is, for
#
#  l_approx(x) = max{ l(x_init) + l'(x_init) * (x - x_init)),
#                     inf_x l(x; (a, b)) }
#
# sets x to minimize
#
#  l_approx(x) + norm(x - x_init)^2 / (2 * stepsize).
#
# Returns the minimizing x.
function TruncatedUpdatePoisson(a::Vector{Float64}, b::Int64,
                                x_init::Vector{Float64}, stepsize::Float64)
  # Replace this code to return the correct update
  return x_init;
end

# x = SGUpdatePoisson(a::Vector{Float64}, b::Int64,
#                     x_init::Vector{Float64}, stepsize::Float64)
#
# Computes the stochastic gradient update for Poisson regression.
function SGUpdatePoisson(a::Vector{Float64}, b::Int64,
                         x_init::Vector{Float64}, stepsize::Float64)
  # Replace this code to return the correct update
  return x_init;
end

max_stars_repo_path: Machine Learning Summer School 2019 (London, UK)/tutorials/optimization/poissonreg.jl


text:
#*************************************************************************#
#
# Preprocessing data
#
#*************************************************************************#
function _compute_stats(mdl::StaticModel,
                        yterms::Vector{Vector{Int}},
                        uterms::Vector{Vector{Int}},
                        iddata::IdData)
    # Get std and mean
    y_std = std(iddata.y, 2)
    u_std = std(iddata.u, 2)
    y_mean = mean(iddata.y, 2)
    u_mean = mean(iddata.u, 2)

    # Compute mean μ and standard deviation σ
    # for each element of the model input x
    ninputs = length_input(mdl)
    k = 1
    σ_x = Vector{Float64}(ninputs)
    μ_x  = Vector{Float64}(ninputs)
    for (i, ynterms) in Enumerate(yterms)
        for yn in ynterms
            σ_x[k] = y_std[i]
            μ_x[k] = y_mean[i]
            k += 1
        end
    end

    for (i, unterms) in Enumerate(uterms)
        for un in unterms
            σ_x[k] = u_std[i]
            μ_x[k] = u_mean[i]
            k += 1
        end
    end

    # Compute mean μ and standard deviation σ
    # for each element of the model input z
    σ_z = vec(y_std)
    μ_z = vec(y_mean)

    return σ_x, μ_x, σ_z, μ_z
end

"""
    learn_offset(mdl, yterms, uterms, iddata) -> new_mdl

Learn offset from data and generate new model with the given offset
embbeded.
"""
function learn_offset(mdl::StaticModel,
                      yterms::Vector{Vector{Int}},
                      uterms::Vector{Vector{Int}},
                      iddata::IdData)

    σ_x, μ_x, σ_z, μ_z = _compute_stats(mdl, yterms, uterms, iddata)

    # Build New Series Model
    new_model = SeriesModel([AffineMap(ones(σ_x), -μ_x), mdl, AffineMap(ones(σ_z), μ_z)])

    return new_model
end

"""
    learn_normalization(mdl, yterms, uterms, iddata[, nσ]) -> new_mdl

Learn offset and scale from data and generate new model with the given offset
and scale embbeded. Where `nσ` is the number of standard deviations the data
will be internally divided by (default is 3).
"""
function learn_normalization{R<:Real}(mdl::StaticModel,
                                      yterms::Vector{Vector{Int}},
                                      uterms::Vector{Vector{Int}},
                                      iddata::IdData;
                                      nσ::R=3.0)

    σ_x, μ_x, σ_z, μ_z = _compute_stats(mdl, yterms, uterms, iddata)

    # Replace σ_x, σ_z by ones when needed
    σ_x = [σ <= 1e-10 ? 1.0 : σ for σ in σ_x]
    σ_z = [σ <= 1e-10 ? 1.0 : σ for σ in σ_z]

    # Build New Series Model
    new_model = SeriesModel([AffineMap(1./(nσ*σ_x) , -μ_x./(nσ*σ_x)), mdl,
                             AffineMap((nσ*σ_z), μ_z)])

    return new_model
end


text:
# @Base.kwdef struct Parameters{T<:AbstractFloat} <: AbstractParameters
#     CSL   :: T  = 0.1   # Surface layer fraction
#     Cτ    :: T  = 0.4   # Von Karman constant
#     CNL   :: T  = 6.33  # Non-local flux proportionality constant
#
#     Cstab :: T  = 2.0   # Stable buoyancy flux parameter for wind-driven turbulence
#     Cunst :: T  = 6.4   # Unstable buoyancy flux parameter for wind-driven turbulence
#
#        Cn :: T  = 1.0   # Exponent for effect of stable buoyancy forcing on wind mixing
#     Cmτ_U :: T  = 0.25  # Exponent for effect of unstable buoyancy forcing on wind mixing of U
#     Cmτ_T :: T  = 0.5   # Exponent for effect of unstable buoyancy forcing on wind mixing of T
#     Cmb_U :: T  = 1/3   # Exponent for the effect of wind on convective mixing of U
#     Cmb_T :: T  = 1/3   # Exponent for effect of wind on convective mixing of T
#
#     Cd_U  :: T  = 0.5   # Wind mixing regime threshold for momentum
#     Cd_T  :: T  = 2.5   # Wind mixing regime threshold for tracers
#
#     Cb_U  :: T  = 0.599 # Buoyancy flux parameter for convective turbulence
#     Cb_T  :: T  = 1.36  # Buoyancy flux parameter for convective turbulence
#     Cτb_U :: T  = (Cτ / Cb_U)^(1/Cmb_U) * (1 + Cunst*Cd_U)^(Cmτ_U/Cmb_U) - Cd_U  # Wind stress parameter for convective turbulence
#     Cτb_T :: T  = (Cτ / Cb_T)^(1/Cmb_T) * (1 + Cunst*Cd_T)^(Cmτ_T/Cmb_T) - Cd_T  # Wind stress parameter for convective turbulence
#
#     CRi   :: T  = 0.3   # Critical bulk Richardson number
#     CKE   :: T  = 4.32  # Unresolved turbulence parameter
#     CKE₀  :: T  = 1e-11 # Minimum unresolved turbulence kinetic energy
#
#     KU₀   :: T  = 1e-6  # Interior viscosity for velocity
#     KT₀   :: T  = 1e-7  # Interior diffusivity for temperature
#     KS₀   :: T  = 1e-9  # Interior diffusivity for salinity
# end

"""
Adapted from
https://github.com/sandreza/OceanConvectionUQSupplementaryMaterials/blob/master/src/ForwardMap/fm.jl
updated for latest version of OceanTurb
"""

"""
closure_kpp_full_evolution(parameters, N, Δt, les::LESbraryData; subsample = 1, grid = 1)

Constructs forward map. Assumes initial conditions and boundary conditions are taken from les data.

# Arguments
- `N`: number of gridpoints to output to
- `Δt`: time step size in seconds
- `les`: les data of the LESbraryData type

# Keyword Arguments
- `subsample`: indices to subsample in time,
- `grid`: in case one wants to save the model grid

# Output
- The forward map. A function that takes parameters and outputs temperature profiles
-   `𝑪`: parameters in KPP, assumes that
    𝑪[1]: Surface Layer Fraction
    𝑪[2]: Nonlocal Flux Amplitude
    𝑪[3]: Diffusivity Amplitude
    𝑪[4]: Shear Constant
"""
function closure_kpp_full_evolution(parameters, U⁰, V⁰, T⁰, t, les; subsample = 1, grid = 1)

     # set parameters
     # parameters = KPP.Parameters( CSL = 𝑪[1], CNL = 𝑪[2], Cb_T = 𝑪[3], CKE = 𝑪[4])

     # assume constant interval between time steps
     Δt = t[2] - t[1]

     # number of gridpoints
     N = length(T⁰)

     # Build the model with a Backward Euler timestepper
     constants = Constants(Float64; α = les.α , β = les.β, f=les.f⁰, g=les.g)
     model = KPP.Model(grid = UniformGrid(N, les.L), stepper=:BackwardEuler, constants = constants, parameters = parameters)

     # Get grid if necessary
     if grid != 1
         zp = collect(model.grid.zc)
         @. grid  = zp
     end

     # Set boundary conditions
     model.bcs.U.top    = FluxBoundaryCondition(les.u_top)
     model.bcs.U.bottom = GradientBoundaryCondition(les.u_bottom)
     model.bcs.T.top    = FluxBoundaryCondition(les.θ_top)
     model.bcs.T.bottom = GradientBoundaryCondition(les.θ_bottom) # may need to fix

    # define the closure
    function evolve()
        # set equal to initial condition of parameterization
        model.solution.U[1:N] = copy(U⁰)
        model.solution.V[1:N] = copy(V⁰)
        model.solution.T[1:N] = copy(T⁰)

        # set aside memory
        if subsample != 1
            time_index = subsample
        else
            time_index = 1:length(t)
        end
        Nt = length(t[time_index])

        U_evolution = zeros(N, Nt)
        V_evolution = zeros(N, Nt)
        T_evolution = zeros(N, Nt)

        # loop the model
        ti = collect(time_index)
        for i in 1:Nt
            time = t[ti[i]]
            run_until!(model, Δt, time)
            @. U_evolution[:,i] = model.solution.U[1:N]
            @. V_evolution[:,i] = model.solution.V[1:N]
            @. T_evolution[:,i] = model.solution.T[1:N]
        end
        return (U_evolution, V_evolution, T_evolution)
    end
    return evolve
end


text:
using PlyIO
using StaticArrays
using Test

@testset "PlyIO" begin

@testset "types" begin
    ply = Ply()
    push!(ply, PlyComment("PlyComment"))
    elt = PlyElement("A", ArrayProperty("x", UInt8[1,2,3]),
                          ArrayProperty("y", Float32[1.1,2.2,3.3]))
    push!(ply, elt)
    @test sprint(show, ply) == """
          Ply with header:
           ply
           format ascii 1.0
           comment PlyComment
           element A 3
           property uint8 x
           property float32 y
           end_header"""
    @test sprint(show, elt) == "PlyElement \"A\" of length 3 with properties [\"x\", \"y\"]"
end

@testset "simple" begin
    ply = Ply()
    push!(ply, PlyComment("PlyComment about A"))

    push!(ply, PlyElement("A",
                          ArrayProperty("x", UInt8[1,2,3]),
                          ArrayProperty("y", Float32[1.1,2.2,3.3]),
                          ListProperty("a_list", Vector{Int64}[[0,1], [2,3,4], [5]])))
    push!(ply, PlyComment("PlyComment about B"))
    push!(ply, PlyComment("PlyObjInfo", obj_info=true))
    push!(ply, PlyComment("PlyComment about B 2"))
    push!(ply, PlyElement("B",
                          ArrayProperty("r", Int16[-1,1]),
                          ArrayProperty("g", Int16[1,1])))
    push!(ply, PlyComment("Final comment"))

    buf = IOBuffer()
    save_ply(ply, buf, ascii=true)
    str = String(take!(buf))
    open("simple_test_tmp.ply", "w") do fid
        write(fid, str)
    end
    @test str ==
    """
    ply
    format ascii 1.0
    comment PlyComment about A
    element A 3
    property uint8 x
    property float32 y
    property list int32 int64 a_list
    comment PlyComment about B
    obj_info PlyObjInfo
    comment PlyComment about B 2
    element B 2
    property int16 r
    property int16 g
    comment Final comment
    end_header
    1	1.1	2 0 1
    2	2.2	3 2 3 4
    3	3.3	1 5
    -1	1
    1	1
    """
end

@testset "empty elements" begin
    ply = Ply()

    push!(ply, PlyElement("A", ArrayProperty("x", UInt8[])))
    push!(ply, PlyElement("B", ListProperty("a_list", Vector{Int64}[[], [], []])))

    buf = IOBuffer()
    save_ply(ply, buf, ascii=true)
    str = String(take!(buf))
    open("empty_test_tmp.ply", "w") do fid
        write(fid, str)
    end
    @test str ==
    """
    ply
    format ascii 1.0
    element A 0
    property uint8 x
    element B 3
    property list int32 int64 a_list
    end_header
    0
    0
    0
    """
end

@testset "roundtrip" begin
    @testset "ascii=$test_ascii" for test_ascii in [false, true]
        ply = Ply()

        push!(ply, PlyComment("A comment"))
        push!(ply, PlyComment("Blah blah"))
        push!(ply, PlyComment("x=10", obj_info=true))

        nverts = 10

        x = collect(Float64, 1:nverts)
        y = collect(Int16, 1:nverts)
        push!(ply, PlyElement("vertex", ArrayProperty("x", x),
                                        ArrayProperty("y", y)))

        # Some triangular faces
        vertex_index = ListProperty("vertex_index", Int32, Int32)
        for i=1:nverts
            push!(vertex_index, rand(0:nverts-1,3))
        end
        push!(ply, PlyElement("face", vertex_index))

        save_ply(ply, "roundtrip_test_tmp.ply", ascii=test_ascii)

        newply = load_ply("roundtrip_test_tmp.ply")

        # TODO: Need a better way to access the data arrays than this.
        @test newply["vertex"]["x"] == x
        @test newply["vertex"]["y"] == y
        @test newply["face"]["vertex_index"] == vertex_index

        @test newply.comments == [PlyComment("A comment",false,1),
                                  PlyComment("Blah blah",false,1),
                                  PlyComment("x=10",true,1)]
    end

    @testset "proptype=$proptype" for proptype in [Int8, Int16, Int32, Int64,
                                                   UInt8, UInt16, UInt32, UInt64,
                                                   Float32, Float64]
        ply = Ply()
        arrayprop = fill(proptype(42), 1)
        listprop = ListProperty("listprop", proptype<:Integer ? proptype : Int32, proptype)
        push!(listprop, collect(proptype, 1:10))
        push!(ply, PlyElement("test", ArrayProperty("arrayprop", arrayprop), listprop))
        io = IOBuffer()
        save_ply(ply, io)
        seek(io, 0)
        newply = load_ply(io)
        @test length(newply) == 1
        @test typeof(newply["test"]["arrayprop"][1]) == proptype
        @test newply["test"]["arrayprop"] == arrayprop
        @test typeof(newply["test"]["listprop"][1]) == Vector{proptype}
        @test newply["test"]["listprop"] == listprop
    end

    @testset "Batched writes for homogenous properties" begin
        ply = Ply()
        nverts = 1000
        x, y, z = [rand(nverts) for _ = 1:3]
        push!(ply, PlyElement("vertex",
                              ArrayProperty("x", x),
                              ArrayProperty("y", y),
                              ArrayProperty("z", z)))
        io = IOBuffer()
        save_ply(ply, io)
        seek(io, 0)
        newply = load_ply(io)
        @test length(newply) == 1
        @test newply["vertex"]["x"] == x
        @test newply["vertex"]["y"] == y
        @test newply["vertex"]["z"] == z
    end
end

@testset "SVector properties" begin
    ply = Ply()
    push!(ply, PlyElement("A",
                          ArrayProperty(["x","y"], SVector{2,Float64}[SVector(1,2), SVector(3,4)])
                         ))
    push!(ply, PlyElement("B",
                          ArrayProperty(["r","g","b"], SVector{3,UInt8}[SVector(1,2,3)])
                         ))
    buf = IOBuffer()
    save_ply(ply, buf, ascii=true)
    str = String(take!(buf))
    open("SVector_properties_test_tmp.ply", "w") do fid
        write(fid, str)
    end
    @test str ==
    """
    ply
    format ascii 1.0
    element A 2
    property float64 x
    property float64 y
    element B 1
    property uint8 r
    property uint8 g
    property uint8 b
    end_header
    1.0	2.0
    3.0	4.0
    1	2	3
    """
end

@testset "Malformed ply headers" begin
    @test_throws ErrorException load_ply(IOBuffer("asdf"))

    @test_throws ErrorException load_ply(IOBuffer("ply"))

    @test_throws ErrorException load_ply(IOBuffer("""
                                                  ply
                                                  format ascii 2.0
                                                  """))

    @test_throws ErrorException load_ply(IOBuffer("""
                                                  ply
                                                  format 1.0
                                                  end_header"""))

    @test_throws ErrorException load_ply(IOBuffer("""
                                                  ply
                                                  format ascii 1.0
                                                  asdf
                                                  end_header"""))

    @test_throws ErrorException load_ply(IOBuffer("""
                                                  ply
                                                  format ascii 1.0
                                                  element el
                                                  end_header"""))

    @test_throws ErrorException load_ply(IOBuffer("""
                                                  ply
                                                  format ascii 1.0
                                                  property float x
                                                  end_header"""))

    @test_throws ErrorException load_ply(IOBuffer("""
                                                  ply
                                                  format ascii 1.0
                                                  element el 0
                                                  property
                                                  end_header"""))

    @test_throws ErrorException load_ply(IOBuffer("""
                                                  ply
                                                  format ascii 1.0
                                                  element el 0
                                                  property list
                                                  end_header"""))
end

end # @testset PlyIO


text:
module MemoizedMethods
using MacroTools: isexpr, combinearg, combinedef, namify, splitarg, splitdef, @capture
export @memoize, memories, forget!

const salt = :__UYTBOOUDVEWICBEWNMDO__ #Symbol(:__, join(rand('A':'Z', 20)), :__)

# which($sig) becomes available in Julia 1.6, so here's a workaround
function _which(tt, world=typemax(UInt))
    meth = ccall(:jl_gf_invoke_lookup, Any, (Any, UInt), tt, world)
    if meth !== nothing
        if meth isa Method
            return meth::Method
        else
            meth = meth.func
            return meth::Method
        end
    end
end

# get_world_counter was never exported, so here's a workaround
_get_world_counter() = ccall(:jl_get_world_counter, UInt, ())

# hasproperty becomes available in Julia 1.2
_hasproperty(x, s) = s in propertynames(x)

"""
    @memoize [cache] declaration

Transform any method declaration `declaration` so that calls to the original method are cached by their arguments. When an argument is unnamed, its type is treated as the argument instead. If the function is defined in local scope, each closure will have a separate cache.

`cache` should be an expression which evaluates to a dictionary-like type that supports `get!` and `empty!`, and may depend on the local variables `__Key__` and `__Value__`, which evaluate to syntactically-determined bounds on the required key and value types the cache must support.

If the given cache contains values, it is assumed that they will agree with the values the method returns. Specializing a method will not empty the cache, but overwriting a method will. You can look up or empty the cache with the functions `memories` and `forget!`.

See also: [`memories`](@ref), [`forget!`](@ref)
"""
macro memoize(args...)
    if length(args) == 1
        cache_constructor = :(IdDict())
        ex = args[1]
    elseif length(args) == 2
        (cache_constructor, ex) = args
    else
        error("Memoize accepts at most two arguments")
    end

    def = try
        splitdef(ex)
    catch
        error("@memoize must be applied to a method definition")
    end

    function split(arg, iskwarg=false)
        arg_name, arg_type, slurp, default = splitarg(arg)
        trait = arg_name === nothing
        trait && (arg_name = gensym())
        vararg = namify(arg_type) === :Vararg
        return (
            arg_name = arg_name,
            arg_type = arg_type,
            arg_value = arg_name,
            slurp = slurp,
            vararg = vararg,
            default = default,
            trait = trait,
            iskwarg = iskwarg)
    end

    combine(arg) = combinearg(arg.arg_name, arg.arg_type, arg.slurp, arg.default)

    pass(arg) = (arg.slurp || arg.vararg) ? Expr(:..., arg.arg_name) :
        arg.iskwarg ? Expr(:kw, arg.arg_name, arg.arg_name) : arg.arg_name

    dispatch(arg) = arg.slurp ? :(Vararg{$(arg.arg_type)}) : arg.arg_type

    key(arg) = arg.trait ? arg.arg_type : arg.arg_name

    key_type(arg) = arg.trait ? DataType :
        arg.vararg ? :(Tuple{$(arg.arg_type)}) : arg.arg_type

    args = split.(def[:args])
    kwargs = split.(def[:kwargs], true)
    def[:args] = combine.(args)
    def[:kwargs] = combine.(kwargs)
    @gensym inferrable
    inferrable_def = deepcopy(def)
    inferrable_args = copy(args)
    inferrable_kwargs = copy(kwargs)
    pop!(inferrable_def, :params, nothing)
    @gensym result

    # If this is a method of a callable type or object, the definition returns nothing.
    # Thus, we must construct the type of the method on our own.
    # We also need to pass the object to the inferrable function
    if haskey(def, :name)
        if haskey(def, :params) # Callable type
            typ = :($(def[:name]){$(pop!(def, :params)...)})
            inferrable_args = [split(:(::Type{$typ})), inferrable_args...]
            def[:name] = combine(inferrable_args[1])
            head = :(Type{$typ})
            sig = :(Tuple{$head, $(dispatch.(args)...)} where {$(def[:whereparams]...)})
        elseif @capture(def[:name], obj_::obj_type_ | ::obj_type_) # Callable object
            inferrable_args = [split(def[:name]), inferrable_args...]
            def[:name] = combine(inferrable_args[1])
            head = obj_type
            sig = :(Tuple{$head, $(dispatch.(args)...)} where {$(def[:whereparams]...)})
        else # Normal syntax style.
            name = def[:name]
            #Named local function definitions are sometimes evaluated
            #asynchronously, and only use Symbols as the function name
            #in both the local and global cases, its safe to name the
            #inference function after the function name, since they
            #dispatch identically. Naming them after the same function
            #name ensures that the symbol is defined when it needs to be.
            if def[:name] isa Symbol
                inferrable = Symbol(def[:name], :_inferrable, salt)
            end

            sig = :(Base.signature_type($(def[:name]), (Tuple{$(dispatch.(args)...)} where {$(def[:whereparams]...)})))
        end
    else # Anonymous function
        head = :(typeof($result))
    end
    tail = :(Tuple{$(dispatch.(args)...)} where {$(def[:whereparams]...)})

    inferrable_def[:name] = inferrable
    inferrable_def[:args] = combine.(inferrable_args)

    cache = gensym(Symbol(:cache, salt))
    bank = Symbol(:bank, salt)

    @gensym getter
    def[:body] = quote
        begin
            $getter = () -> $inferrable($(pass.(inferrable_args)...); $(pass.(inferrable_kwargs)...))
            $(get!)($getter, $cache[2], ($(map(key, [inferrable_args; inferrable_kwargs])...),))
        end::$(Core.Compiler.widenconst)($(Core.Compiler.return_type)($getter, $(Tuple{})))
    end

    scope = gensym()

    quote
        # The `local` qualifier will make this performant even in the global scope.
        $(esc(quote
            local $cache = begin
                local __Key__ = (Tuple{$(map(key_type, [inferrable_args; inferrable_kwargs])...)} where {$(def[:whereparams]...)})
                local __Value__ = ($(get(def, :rtype, Any)) where {$(def[:whereparams]...)})
                ($tail, $cache_constructor)
            end
        end))

        $(if @isdefined sig
            quote
                local world = _get_world_counter()
            end
        end)

        $(esc(combinedef(inferrable_def)))
        local $(esc(result)) = Base.@__doc__($(esc(combinedef(def))))

        $(if @isdefined sig
            quote
                if !_hasproperty($(esc(result)), $(Expr(:quote, cache)))
                    local sig = $(esc(sig))
                    # If overwriting a method, empty the old cache.
                    # Notice that methods are hashed by their stored signature
                    local meth = _which(sig, world)
                    if meth != nothing && meth.sig == sig && isdefined(meth.module, $(Expr(:quote, bank)))
                        empty!(pop!(meth.module.$bank, meth.sig, (nothing, []))[2])
                    end

                    if !isdefined($__module__, $(Expr(:quote, bank)))
                        global $(esc(bank)) = IdDict()
                    end
                    # Store the cache so that it can be emptied later
                    local meth = _which(sig)
                    $(esc(bank))[meth.sig] = $(esc(cache))
                end
            end
        end)

        $(esc(result))
    end
end

"""
    memories(f, types::Type)

If the method `which(f, types)`, is memoized, return the cache in the
scope of `f`. Otherwise, return `nothing` or an overwritten cache.

See also: [`@memoize`](@ref), [`forget!`](@ref)
"""
function memories(f, types)
    m = which(f, types)
    types = Tuple{m.sig.parameters[2:end]...}
    for name in propertynames(f) #if f is a closure, we walk its fields
        if first(string(name), length(string("##cache", salt))) == string("##cache", salt)
            cache = getproperty(f, name)
            if cache isa Core.Box
                cache = cache.contents
            end
            if (cache[1] == types)
                return cache[2]
            end
        end
    end
    return memories(m) #otherwise, a method would suffice
end

"""
    memories(m::Method)

If m is a memoized method defined at global scope, return its cache.
Otherwise, return `nothing` or an overwritten cache.

See also: [`@memoize`](@ref), [`forget!`](@ref)
"""
function memories(m::Method)
    if isdefined(m.module, Symbol(:bank, salt))
        return get(getproperty(m.module, Symbol(:bank, salt)), m.sig, (nothing, nothing))[2]
    end
    return nothing
end

"""
    forget!(f, types::Type)

If the method `which(f, types)`, is memoized, `empty!` the cache in the
scope of `f`.

See also: [`@memoize`](@ref), [`memories`](@ref)
"""
function forget!(f, types::Type)
    c = memories(f, types)
    if c !== nothing
        return empty!(c)
    else
        return nothing
    end
end

"""
    forget!(m::Method)

If the method `m`, is memoized, `empty!` its cache. To clear all memoized
caches of a global method `f`, you can `map(forget!, methods(f))`.

See also: [`@memoize`](@ref), [`memories`](@ref)
"""
function forget!(m::Method)
    c = memories(m)
    if c !== nothing
        return empty!(c)
    else
        return nothing
    end
end

end


text:
import IncrementalInference: getSample
import Base: convert, show
# import DistributedFactorGraphs: getManifold

using .Images

export ScanMatcherPose2, PackedScanMatcherPose2

"""
$TYPEDEF

This is but one incarnation for how radar alignment factor could work, treat it as a starting point.

Notes
- Stanard `cvt` argument is lambda function to convert incoming images to user convention of image axes,
  - default `cvt` flips image rows so that Pose2 xy-axes corresponds to img[x,y] -- i.e. rows down and across from top left corner.
- Use rescale to resize the incoming images for lower resolution (faster) correlations
- Both images passed to the construct must have the same type some matrix of type `T`.

Example
-------
```julia
arp2 = ScanMatcherPose2(img1, img2, 2) # e.g. 2 meters/pixel
```

See also: [`overlayScanMatcher`](@ref)
"""
struct ScanMatcherPose2{T} <: IIF.AbstractManifoldMinimize
  """ reference image for scan matching. """
  im1::Matrix{T}
  """ test image to scan match against the reference image. """
  im2::Matrix{T}
  """ Common grid scale for both images -- i.e. units/pixel.
  Constructor uses two arguments `gridlength`*`rescale=1`=`gridscale`.
  Arg 0 < `rescale` ≤ 1 is also used to rescale the images to lower resolution for speed. """
  gridscale::Float64

  # replace inner constructor with transform on image
  ScanMatcherPose2{T}( im1::AbstractMatrix{T},
                      im2::AbstractMatrix{T},
                      gridlength::Real,
                      rescale::Real=1,
                      cvt = (im)->reverse(Images.imresize(im,trunc.(Int, rescale.*size(im))),dims=1)
                    ) where {T} = new{T}( cvt(im1),
                                          cvt(im2),
                                          rescale*gridlength )
end

ScanMatcherPose2(im1::AbstractMatrix{T}, im2::AbstractMatrix{T}, w...) where T = ScanMatcherPose2{T}(im1,im2, w...)

getManifold(::IIF.InstanceType{<:ScanMatcherPose2}) = getManifold(Pose2Pose2)

getSample( cf::CalcFactor{<:ScanMatcherPose2} ) = nothing

function (cf::CalcFactor{<:ScanMatcherPose2})(X, p, q)
  M = getManifold(Pose2)
  arp = cf.factor
  # for groups
  tf = Manifolds.compose(M, inv(M, p), q)
  Mr = M.manifold[2]
  r0 = identity_element(Mr)
  r = vee(Mr, r0, log(Mr, r0, tf.parts[2]))[1]
  return evaluateTransform(arp.im1, arp.im2, tf.parts[1], r, arp.gridscale)
end

function Base.show(io::IO, arp::ScanMatcherPose2{T}) where {T}
  printstyled(io, "ScanMatcherPose2{", bold=true, color=:blue)
  println(io)
  printstyled(io, "    T = ", color=:magenta)
  println(io, T)
  printstyled(io, " }", color=:blue, bold=true)
  println(io)
  println(io, "  size(.im1):      ", size(arp.im1))
  println(io, "    min/max:       ", round(minimum(arp.im1),digits=3), "/", round(maximum(arp.im1),digits=3))
  println(io, "  size(.im2):      ", size(arp.im2))
  println(io, "    min/max:       ", round(minimum(arp.im2),digits=3), "/", round(maximum(arp.im2),digits=3))
  println(io, "  gridscale:       ", arp.gridscale)
  nothing
end

Base.show(io::IO, ::MIME"text/plain", arp::ScanMatcherPose2) = show(io, arp)
Base.show(io::IO, ::MIME"application/juno.inline", arp::ScanMatcherPose2) = show(io, arp)

"""
    $SIGNATURES

Overlay the two images from `AlignRadarPose2` with the first (red) fixed and transform the second image (blue) according to `tf`.

Notes:
- `tf` is a Manifolds.jl type `::ProductRepr` (or newer `::ArrayPartition`) to represent a `SpecialEuclidean(2)` manifold point.
"""
function overlayScanMatcher(sm::ScanMatcherPose2,
                            trans::AbstractVector{<:Real}=Float64[0;0],
                            rot::Real=0.0;
                            score=Ref(0.0),
                            showscore::Bool=true  )
  #
  im2_ = transformImage_SE2(sm.im2, trans, rot, sm.gridscale)

  # get matching padded views
  im1_, im2__ = paddedviews(0.0, sm.im1, im2_)

  im1__ = RGBA.(im1_, 0, 0, 0.5)
  im2___ = RGBA.(0, 0, im2__, 0.5)

  if showscore
    score[] = evaluateTransform(sm.im1, sm.im2, trans, rot, sm.gridscale)
    @info "overlayScanMatcher score" score[]
  end

  # im2__ = RGBA.(im2_, 0, 0, 0.5)
  im1__ .+ im2___
end

#

function overlayScanMatcher(sm::ScanMatcherPose2,
                            tf = Manifolds.identity_element(SpecialEuclidean(2));
                            kw... )
  #
  M = SpecialOrthogonal(2)
  e0 = identity_element(M)

  rot = vee( M, e0, log(M, e0, tf.parts[2]) )[1]

  overlayScanMatcher(sm, tf.parts[1], rot; kw...)
end

## =========================================================================================
## Factor serialization below
## =========================================================================================

struct PackedScanMatcherPose2 <: PackedInferenceType
  im1::Vector{Vector{Float64}}
  im2::Vector{Vector{Float64}}
  gridscale::Float64
end

function convert(::Type{<:PackedScanMatcherPose2}, arp2::ScanMatcherPose2)
  TensorCast.@cast pim1[row][col] := arp2.im1[row,col]
  TensorCast.@cast pim1[row] := collect(pim1[row])
  TensorCast.@cast pim2[row][col] := arp2.im2[row,col]
  TensorCast.@cast pim2[row] := collect(pim2[row])
  PackedScanMatcherPose2(
    pim1,
    pim2,
    arp2.gridscale )
end

function convert(::Type{<:ScanMatcherPose2}, parp2::PackedScanMatcherPose2)
  TensorCast.@cast im1[row,col] := parp2.im1[row][col]
  TensorCast.@cast im2[row,col] := parp2.im2[row][col]
  ScanMatcherPose2(
    collect(im1),
    collect(im2),
    parp2.gridscale )
end

#


text:
import Effects: effects
import Effects:_reference_grid
import Effects:typify
import Effects.typify
import Effects:_symequal
import StatsModels.collect_matrix_terms
import Base.getproperty
using Effects

"""
effects(design::AbstractDict, model::UnfoldModel;typical=mean)

Calculates marginal effects for all term-combinations in `design`.

 Implementation based on Effects Package; likely could repackage in UnfoldEffects; somebody wants to do it? This would make it easier to cross-maintain it to changes/bugfixes in the Effects.jl Package
 `design` is a Dictionary containing those predictors (als keys) with levels (as values), that you want to evaluate. The `typical` refers to the value, that other predictors not in the Dictionary should take on.

 # Example
 ```julia-repl
 julia> f = @formula 0 ~ categoricalA + continuousA + continuousB
 julia> uf = fit(UnfoldModel,(Any=>(f,times)),data,events)
 julia> d = Dict(:categoricalA=>["levelA","levelB"],:continuousB=>[-2,0,2])
 julia> effects(d,uf)
```
 will result in 6 predicted values: A/-2, A/0, A/2, B/-2, B/0, B/2.
"""

function effects(design::AbstractDict, model::UnfoldModel;typical=mean)
    reference_grid = _reference_grid(design)
    form = Unfold.formula(model) # get formula

    # replace non-specified fields with "constants"
    m = modelmatrix(model,false) # get the modelmatrix without timeexpansion
    form_typical = _typify(reference_grid,form,m,typical)

    eff = yhat(model,form_typical,reference_grid)

    # because coefficients are 2D/3D arry, we have to cast it correctly to one big dataframe
    if isa(eff,Tuple)
        # TimeContinuous Model, we also get back other things like times & fromWhereToWhere a BasisFunction goes
        if typeof(form) <: FormulaTerm
            form = [form]
        end
        # figure out the basisname
        bnames = [[form[ix].rhs.basisfunction.name] for ix in 1:length(form)]
        # repeat it as much as necessary
        bnames = repeat.(bnames,[e.stop+e.step-1 for e in eff[1]])

        result = DataFrame(cast_referenceGrid(reference_grid,eff[3],eff[2] ;basisname=vcat(bnames...)))
        select!(result,Not(:latency)) # remove the latency column if it was added
    else
        # normal mass univariate model
        result = DataFrame(cast_referenceGrid(reference_grid,eff,times(model)[1] ))
    end

return result
end

 Effects.typify(reference_grid,form::Matrix,X;kwargs...) = typify.(Ref(reference_grid),form,Ref(X);kwargs...)

 # cast single form to a vector
_typify(reference_grid,form::FormulaTerm{<:InterceptTerm,<:Unfold.TimeExpandedTerm},m,typical) = _typify(reference_grid,[form],[m],typical)

function _typify(reference_grid, form::AbstractArray{<:FormulaTerm{<:Union{<:InterceptTerm,<:Unfold.TimeExpandedTerm}}},m,typical)

    form_typical = Array{Any}(undef,1, length(form))
    for f = 1:length(form)

        # strip of basisfunction and put it on afterwards again
        tmpf = deepcopy(form[f])

        # create a Formula without Timeexpansion
        tmpf = FormulaTerm(tmpf.lhs,tmpf.rhs.term)

        # typify that
        tmpf = typify(reference_grid, tmpf, m[f]; typical=typical)

        # regenerate TimeExpansion
        tmpf = Unfold.TimeExpandedTerm(tmpf,form[f].rhs.basisfunction;eventfields=form[f].rhs.eventfields)
        form_typical[f] = tmpf
    end
    return form_typical
 end
 function _typify(reference_grid,form::FormulaTerm,m,typical)
    return [typify(reference_grid, form, m; typical=typical)]

 end

function cast_referenceGrid(r,eff,times;basisname=nothing)
    nchan = size(eff, 2) # correct
    neff = size(r,1) # how many effects requested
    neffCol = size(r,2) # how many predictors
    ncols = size(eff,1) ÷ neff # typically ntimes

    # replicate
    # for each predictor in r (reference grid), we need this at the bottom

    if isnothing(basisname)
        nbases = 1
    else
        nbases = length(unique(basisname))
    end
    coefs_rep = Array{Array}(undef,nbases,neffCol)

    for k = 1:neffCol
        # in case we have only a single basis (e.g. mass univariate), we can directly fill in all values
        ixList = []
        if isnothing(basisname)
            ix = ones(ncols) .== 1.
            append!(ixList,[ix])
        else
            #in case of multiple bases, we have to do it iteratively, because the bases can be different length
            for b = unique(basisname)
                ix = basisname[1:neff:end].==b
                append!(ixList,[ix])
            end
        end
        for i_ix = 1:length(ixList)

            coefs_rep[i_ix,k] = linearize(permutedims(repeat(r[:,k], outer = [1, nchan, sum(ixList[i_ix])]), [2,3,1]))
        end
    end

    # often the "times" vector
    if length(times)  == neff*ncols
        # in case we have timeexpanded, times is already in long format and doesnt need to be repeated for each coefficient
        colnames_basis_rep = permutedims(repeat(times, 1, nchan, 1), [2 1 3])
    else
        colnames_basis_rep = permutedims(repeat(times, 1, nchan, neff), [2 1 3])
    end

    # for multiple channels
    chan_rep = repeat(1:nchan, 1, ncols, neff)

    # for mass univariate there is no basisname
    if isnothing(basisname)
        basisname = fill(nothing,ncols)
        basisname_rep = permutedims(repeat(basisname,1,nchan,neff),[2,1,3])
    else
        basisname_rep = permutedims(repeat(basisname,1,nchan,1),[2,1,3])
    end

    result = Dict(   :yhat => linearize(eff'),
            :time => linearize(colnames_basis_rep),
            :channel => linearize(chan_rep),
            :basisname =>linearize(basisname_rep))

    for k = 1:neffCol
            push!(result,Symbol(names(r)[k]) =>vcat(vcat(coefs_rep[:,k]...)...))
    end

    return result
end
Effects._trmequal(t1::Unfold.uf_bsplineTerm,t2::AbstractTerm) = _symequal(t1.term,t2)
Effects._trmequal(t1::Unfold.uf_bsplineTerm,t2::Unfold.uf_bsplineTerm) = _symequal(t1.term,t2.term)

Effects._trmequal(t1::AbstractTerm,t2::Unfold.uf_bsplineTerm) = _symequal(t1,t2.term)
Effects._symequal(t1::AbstractTerm,t2::Unfold.uf_bsplineTerm) = _symequal(t1,t2.term)
#Effects._symequal(t1::AbstractTerm,t2::Unfold.TimeExpandedTerm) = _symequal(t1,t2.term)
#function Effects._replace(matrix_term::MatrixTerm{<:Tuple{<:Unfold.TimeExpandedTerm}},typicals::Dict)

#    replaced_term = MatrixTerm((Effects._replace.(matrix_term.terms, Ref(typicals))...,))
#    basisfunctionTerm = getfield(matrix_term,:terms)[1]
#    return TimeExpandedTerm(replaced_term,basisfunctionTerm.basisfunction,basisfunctionTerm.eventfields)

#end

max_stars_repo_name: unfoldtoolbox/jlunfold

max_issues_repo_name: unfoldtoolbox/jlunfold

max_forks_repo_name: unfoldtoolbox/jlunfold


text:
"""
`Maker.DirectoryTarget`

The type created by `directory()`. Fields expected to be accessed publicly
include:

- `name::UTF8String`
- `dependencies::Vector{UTF8String}`
- `description::UTF8String`

"""
type DirectoryTarget <: AbstractTarget
    name::UTF8String
    dependencies::Vector{UTF8String}
    description::UTF8String
    action::Function
    timestamp::DateTime
    funhash::UInt64
    isstale::Bool
end

function isstale(t::DirectoryTarget)
    if !isdir(t.name) || t.isstale
        return true
    end
    ds = dependencies(t)
    length(ds) == 0 ? false : timestamp(t) < maximum([timestamp(d) for d in ds])
end

function execute(t::DirectoryTarget)
    !isdir(t.name) && mkpath(t.name)
    t.isstale = false
    updatecache!(t)
end

timestamp(t::DirectoryTarget) = if isdir(t.name)
    Dates.unix2datetime(0)
else
    Dates.unix2datetime(time())  # this will trigger the target
end

"""
```julia
directory(name::AbstractString, dependencies=[])
```
Define and register targets for Maker.jl.

- `name` refers to the name of the task or target.

- `dependencies` refers to the name (`AbstractString`) or names
  (`Vector{AbstractString}`) of targets that need to be satisfied
  for this target.

Targets are registered globally.

`directory` targets use the name of the path as the name of the target.
No action can be specified. The path `name` is created with `mkpath` if
it doesn't exist. The path is relative to the current working directory.

`directory` targets only have a builtin action that creates the
directory if it is missing. A timestamp is stored when the directory
is created, but the return from `timestamp("directorytask")` returns
the oldest time if the directory exists (it doesn't need to be recreated),
and it returns the current time if the directory doesn't exist (it needs
 to be created).

See also `make`, `file`, and `task`. `directory` registers a
`DirectoryTarget` type.
"""
directory(name::AbstractString, dependencies=UTF8String[]) =
    target(DirectoryTarget, name, () -> nothing, dependencies)


text:
import StatsBase: ProbabilityWeights, pweights

v1 = UncertainValue(UnivariateKDE, rand(4:0.25:6, 1000), bandwidth = 0.01)
v2 = UncertainValue(Normal, 0.8, 0.4)
v3 = UncertainValue([rand() for i = 1:3], [0.3, 0.3, 0.4])
v4 = UncertainValue(Normal, 3.7, 0.8)
uvals = [v1, v2, v3, v4];

# Combining without weights
r1 = combine(uvals)
r2 = combine(uvals, n = 10000)
@test r1 isa AbstractUncertainValue
@test r2 isa AbstractUncertainValue

# Combining with ProbabilityWeights
r1 = combine(uvals, ProbabilityWeights([0.2, 0.1, 0.3, 0.2]))
r2 = combine(uvals, pweights([0.2, 0.1, 0.3, 0.2]))
@test r1 isa AbstractUncertainValue
@test r2 isa AbstractUncertainValue

# Combining with AnalyticWeights
combine(uvals, AnalyticWeights([0.2, 0.1, 0.3, 0.2]))
combine(uvals, aweights([0.2, 0.1, 0.3, 0.2]))
@test r1 isa AbstractUncertainValue
@test r2 isa AbstractUncertainValue

# Combining with FrequencyWeights
combine(uvals, FrequencyWeights([100, 200, 300, 400]))
combine(uvals, fweights([500, 700, 800, 124]))
@test r1 isa AbstractUncertainValue
@test r2 isa AbstractUncertainValue

# Combining with generic Weights
r1 = combine(uvals, Weights([0.2, 0.1, 0.3, 0.2]))
@test r1 isa AbstractUncertainValue


text:
# ------------------------------------
# Tests adapted from Git2Go Library
# ------------------------------------
#=
context("odb test") do
    test_path = joinpath(pwd(), "testrepo")
    repo = create_test_repo(test_path)
    try
        cid, tid = seed_test_repo(repo)
        str = "hello, world!"
        odb = Odb(repo)
        stream = open_wstream(GitBlob, odb, length(str))
        @test iswriteable(stream)
        @test !isreadable(stream)
        write(stream, str)
        close(stream)
        @show eid = Oid("30f51a3fba5274d53522d0f19748456974647b4f")
        @test eid == Oid(stream)
    finally
        close(repo)
        cleanup_dir(test_path)
    end
end
=#


text:
"""
    A collection of analytic tools to evaluate CellML models
"""

function simplify_systems(systems)
    for x in systems
        sys = last(x)
        print("simplifying $(nameof(sys))")
        try
            sys = structural_simplify(sys)
            k = max(1, 50-length(string(nameof(sys))))
            printstyled(repeat(" ", k) * "OK!"; color=:green)
            println()
        catch e
            println(e)
        end

    end
end