text:
import IntertemporalChoiceHeuristics:
    load,
    fit,
    models,
    evaluation_losses,
    average_loss

import Optim: minimizer

inputs, weights = load()

for model in models
    res = fit(model, inputs, weights)
    parameters = minimizer(res)
    m = model(parameters)
    for loss in evaluation_losses
        l = average_loss(m, loss, inputs, weights)
        @printf(
            "%s\t%s\t%s\t%s\t%s\n",
            m,
            loss,
            l,
            Optim.minimum(res),
            norm(parameters),
        )
    end
end

ext: jl


text:
export lanczos

# this finds the OPs and recurrence for a
function lanczos(w,N)
    x = Fun(identity,space(w))

    f1=Fun(1/sqrt(sum(w)),space(x))

    P = Array{Fun}(undef, N + 1)
    β = Array{eltype(w)}(undef, N)
    γ = Array{eltype(w)}(undef, N)

    P[1] = f1

    v = x*P[1]
    β[1] = sum(w*v*P[1])

    v = v - β[1]*P[1]
    γ[1] = sqrt(sum(w*v^2))

    P[2] = v/γ[1]

    for k = 2:N
        v = x*P[k] - γ[k-1]*P[k-1]
        β[k] = sum(w*v*P[k])
        v = v - β[k]*P[k]
        γ[k] = sqrt(sum(w*v^2))
        P[k+1] = v/γ[k]
    end

    P,β,γ
end


text:
# # Type Recipes

# this is the default "type recipe"... just pass the object through
@recipe f(::Type{T}, v::T) where {T} = v

# this should catch unhandled "series recipes" and error with a nice message
@recipe f(::Type{V}, x, y, z) where {V <: Val} =
    error("The backend must not support the series type $V, and there isn't a series recipe defined.")

"""
    _apply_type_recipe(plotattributes, v::T, letter)

Apply the type recipe with signature `(::Type{T}, ::T)`.
"""
function _apply_type_recipe(plotattributes, v, letter)
    _preprocess_axis_args!(plotattributes, letter)
    rdvec = RecipesBase.apply_recipe(plotattributes, typeof(v), v)
    warn_on_recipe_aliases!(plotattributes[:plot_object], plotattributes, :type, typeof(v))
    _postprocess_axis_args!(plotattributes, letter)
    return rdvec[1].args[1]
end

# Handle type recipes when the recipe is defined on the elements.
# This sort of recipe should return a pair of functions... one to convert to number,
# and one to format tick values.
function _apply_type_recipe(plotattributes, v::AbstractArray, letter)
    plt = plotattributes[:plot_object]
    _preprocess_axis_args!(plotattributes, letter)
    # First we try to apply an array type recipe.
    w = RecipesBase.apply_recipe(plotattributes, typeof(v), v)[1].args[1]
    warn_on_recipe_aliases!(plt, plotattributes, :type, typeof(v))
    # If the type did not change try it element-wise
    if typeof(v) == typeof(w)
        isempty(skipmissing(v)) && return Float64[]
        x = first(skipmissing(v))
        args = RecipesBase.apply_recipe(plotattributes, typeof(x), x)[1].args
        warn_on_recipe_aliases!(plt, plotattributes, :type, typeof(x))
        _postprocess_axis_args!(plotattributes, letter)
        if length(args) == 2 && all(arg -> arg isa Function, args)
            numfunc, formatter = args
            return Formatted(map(numfunc, v), formatter)
        else
            return v
        end
    end
    _postprocess_axis_args!(plotattributes, letter)
    return w
end

# special handling for Surface... need to properly unwrap and re-wrap
_apply_type_recipe(
    plotattributes,
    v::Surface{<:AMat{<:Union{AbstractFloat, Integer, AbstractString, Missing}}},
    letter,
) = v
function _apply_type_recipe(plotattributes, v::Surface, letter)
    ret = _apply_type_recipe(plotattributes, v.surf, letter)
    if typeof(ret) <: Formatted
        Formatted(Surface(ret.data), ret.formatter)
    else
        Surface(ret)
    end
end

# don't do anything for datapoints or nothing
_apply_type_recipe(
    plotattributes,
    v::AbstractArray{<:Union{AbstractFloat, Integer, AbstractString, Missing}},
    letter,
) = v
_apply_type_recipe(plotattributes, v::Nothing, letter) = v

# axis args before type recipes should still be mapped to all axes
function _preprocess_axis_args!(plotattributes)
    plt = plotattributes[:plot_object]
    for (k, v) in plotattributes
        if is_axis_attribute(plt, k)
            pop!(plotattributes, k)
            for l in (:x, :y, :z)
                lk = Symbol(l, k)
                haskey(plotattributes, lk) || (plotattributes[lk] = v)
            end
        end
    end
end
function _preprocess_axis_args!(plotattributes, letter)
    plotattributes[:letter] = letter
    _preprocess_axis_args!(plotattributes)
end

# axis args in type recipes should only be applied to the current axis
function _postprocess_axis_args!(plotattributes, letter)
    plt = plotattributes[:plot_object]
    pop!(plotattributes, :letter)
    if letter in (:x, :y, :z)
        for (k, v) in plotattributes
            if is_axis_attribute(plt, k)
                pop!(plotattributes, k)
                lk = Symbol(letter, k)
                haskey(plotattributes, lk) || (plotattributes[lk] = v)
            end
        end
    end
end


text:
#!/usr/bin/env julia

# An example for using Gadfly within GTK

using Gtk
using Cairo
using Compose
using Gadfly
using GtkBuilderAid

type PlotData
  surface::CairoSurface
  plot
end

function draw_backbuffer(area, udata)
  udata.surface = create_similar_surface(area, Gtk.GEnum(Cairo.CONTENT_COLOR_ALPHA))
  Gadfly.draw(
    Compose.CAIROSURFACE(udata.surface, CairoContext(udata.surface)),
    udata.plot)
end

@GtkBuilderAid function_name(canvas_builder) begin

  @guarded Cint(1) function plot_area_draw_cb(area, ctx_ptr, udata)
    ctx = CairoContext(ctx_ptr)
    set_source_surface(ctx, udata.surface)
    paint(ctx)
    Cint(0)
  end

  @guarded Cint(0) function plot_area_configure_event_cb(
      area,
      configure_event,
      udata)
    # redraw the plot
    draw_backbuffer(area, udata)
    Cint(1)
  end

  @guarded function plot_area_realize_cb(
      area,
      udata)
    draw_backbuffer(area, udata)
    nothing
  end

end

@guarded function activate_cb(
    app_ptr,
    userdata)
  app = GObject(app_ptr)
  built = canvas_builder(
    "resources/main.ui",
    PlotData(
      CairoSurface(C_NULL, -1, -1),
      plot(x=collect(1:10), y=rand(10), Geom.LineGeometry)))
  win = GAccessor.object(built, "main_window")
  push!(app, win)
  showall(win)
  nothing
end

app = @GtkApplication("io.github.matt5sean3.GtkBuilderAid.fourth", 0)
signal_connect(activate_cb, app, :activate, Void, (), false, ())

println("Starting Application")
run(app)


text:
export cg,bicgstab,Lanczos

function cg{T1<:Number,T2<:Number}(A::AbstractMatrix{T1},b::Vector{T2})
	n=length(b)
	n1,n2=size(A)
	n == n1 == n2 || throw(DimensionMismatch(""))
	T=promote_type(T1,T2)
	x=zeros(T,n)
	r=b-A*x
	p=r
	rsold=dot(r,r)
	TR=typeof(real(rsold))
	for i=1:n
		Ap=A*p
		alpha=rsold/dot(p,Ap)
		x=x+alpha*p
		r=r-alpha*Ap
		rsnew=dot(r,r)
		if real(sqrt(rsnew))<10eps(TR)
			return x
		end
		p=r+rsnew/rsold*p
		rsold=rsnew
	end
	return x
end

function cg{T1<:Number,T2<:Number,T3<:Number}(A::AbstractMatrix{T1},b::Vector{T2},P::AbstractMatrix{T3})
	n=length(b)
	n1,n2=size(A)
	n3,n4=size(P)
	n == n1 == n2 == n3 == n4 || throw(DimensionMismatch(""))
	T=promote_type(T1,T2,T3)
	x=P\b
	r=b-A*x
	z=P\r
	p=z
	rsold=dot(r,z)
	TR=typeof(real(rsold))
	for i=1:n
		Ap=A*p
		alpha=rsold/dot(p,Ap)
		x=x+alpha*p
		r=r-alpha*Ap
		z=P\r
		rsnew=dot(r,z)
		if real(sqrt(rsnew))<10eps(TR)
			return x
		end
		p=z+rsnew/rsold*p
		rsold=rsnew
	end
	return x
end

function bicgstab{T1<:Number,T2<:Number}(A::AbstractMatrix{T1},b::Vector{T2})
	n=length(b)
	n1,n2=size(A)
	n == n1 == n2 || throw(DimensionMismatch(""))
	T=promote_type(T1,T2)
	x=zeros(T,n)
	r=b-A*x
	rhat = r
	ρold = one(T)
	α =one(T)
	ωold = one(T)
	v=zeros(T,n)
	p=zeros(T,n)
	TR = typeof(real(α))
	for i=1:n
		ρnew = dot(rhat,r)
		β = (ρnew/ρold)*(α/ωold)
		p = r + β*(p-ωold*v)
		v = A*p
		α = ρnew/dot(rhat,v)
		s = r - α*v
		t = A*s
		ωnew = dot(t,s)/dot(t,t)
		upd = α*p + ωnew*s
		x = x + upd
		if norm(upd)<10eps(TR)
			return x
		end
		r = s - ωnew*t
		ρold = ρnew
		ωold = ωnew
	end
	return x
end

function bicgstab{T1<:Number,T2<:Number,T3<:Number}(A::AbstractMatrix{T1},b::Vector{T2},P::AbstractMatrix{T3})
	n=length(b)
	n1,n2=size(A)
	n3,n4=size(P)
	n == n1 == n2 == n3 == n4 || throw(DimensionMismatch(""))
	T=promote_type(T1,T2,T3)
	x=P\b
	r=b-A*x
	rhat = r
	ρold = one(T)
	α =one(T)
	ωold = one(T)
	v=zeros(T,n)
	p=zeros(T,n)
	TR = typeof(real(α))
	for i=1:n
		ρnew = dot(rhat,r)
		β = (ρnew/ρold)*(α/ωold)
		p = r + β*(p-ωold*v)
		y = P\p
		v = A*y
		α = ρnew/dot(rhat,v)
		s = r - α*v
		z = P\s
		t = A*z
		zt = P\t
		ωnew = dot(zt,z)/dot(zt,zt)
		upd = α*y + ωnew*z
		x = x + upd
		if norm(upd)<10eps(TR)
			return x
		end
		r = s - ωnew*t
		ρold = ρnew
		ωold = ωnew
	end
	return x
end

function Lanczos{T<:Number}(A::AbstractMatrix{T})
	n,n1 = size(A)
	n == n1 || throw(DimensionMismatch(""))
	z = Array(T,n)
	vold = zero(T)
	vnew = [one(T),zeros(T,n-1)]
	beta = zeros(T,n)
	alpha = zeros(T,n)
	for j = 1:n-1
		w = A*vnew
		alpha[j] = dot(w,vnew)
		w = w - alpha[j]*vnew - beta[j]*vold
		beta[j+1] = norm(w)
		vold = vnew
		vnew = w/beta[j+1]
	end
	w = A*vnew
	alpha[n] = dot(w,vnew)
	@time steig!(alpha,[beta[2:end],zero(T)],z,10000)
	return sort!(alpha)
 end


text:
type CellSpace
    cells::Array{Cell,2} # cell array
    width::Integer         # cell space width
    height::Integer        # cell space height
    function CellSpace(w::Integer, h::Integer)
        cells = zeros(Cell, w, h)
        new(cells, w, h)
    end
    function CellSpace(cells::Array{Cell,2})
        w = size(cells, 1)
        h = size(cells, 2)
        new(cells, w, h)
    end
end

# 0x11223344 <-> 0xNNEESSWW
@inline function north(val::Cell)
    (val >> BSHIFT_NORTH) & BMASK_P
end
@inline function east(val::Cell)
    (val >> BSHIFT_EAST) & BMASK_P
end
@inline function south(val::Cell)
    (val >> BSHIFT_SOUTH) & BMASK_P
end
@inline function west(val::Cell)
    (val >> BSHIFT_WEST) & BMASK_P
end

@inline function qn(d::SubArray)
    Transition(south(d[QNORTH_IDX]))
end
@inline function qe(d::SubArray)
    Transition(west(d[QEAST_IDX]))
end
@inline function qs(d::SubArray)
    Transition(north(d[QSOUTH_IDX]))
end
@inline function qw(d::SubArray)
    Transition(east(d[QWEST_IDX]))
end

@inline function qn!(d::SubArray, val::Cell)
    d[QNORTH_IDX] = (d[QNORTH_IDX] & ~BMASK_SOUTH)|(val << BSHIFT_SOUTH)
end
@inline function qe!(d::SubArray, val::Cell)
    d[QEAST_IDX] = (d[QEAST_IDX] & ~BMASK_WEST)|(val << BSHIFT_WEST)
end
@inline function qs!(d::SubArray, val::Cell)
    d[QSOUTH_IDX] = (d[QSOUTH_IDX] & ~BMASK_NORTH)|(val << BSHIFT_NORTH)
end
@inline function qw!(d::SubArray, val::Cell)
    d[QWEST_IDX] = (d[QWEST_IDX] & ~BMASK_EAST)|(val << BSHIFT_EAST)
end

function get_target_state(d::SubArray)
    (Transition(d[TARGET_IDX]) << BSHIFT_CELL) |
    (qn(d) << BSHIFT_NORTH) |
    (qe(d) << BSHIFT_EAST) |
    (qs(d) << BSHIFT_SOUTH) |
    (qw(d) << BSHIFT_WEST)
end

function set_target_state(d::SubArray, val::Transition)
    hval = Cell(val >> BSHIFT_CELL)
    lval = Cell(val & BMASK_C)
    d[TARGET_IDX] = hval
    qn!(d, north(lval))
    qe!(d, east(lval))
    qs!(d, south(lval))
    qw!(d, west(lval))
end

function update_checkerboard!(cs::CellSpace, rule::Rule)
    w = cs.width-1
    h = cs.height-1

    @parallel for idx in [vec([(i,j) for i=3:2:w,j=3:2:h]);vec([(i,j) for i=2:2:w,j=2:2:h])]
        update!(cs, rule, idx[1], idx[2])
    end
    @parallel for idx in [vec([(i,j) for i=2:2:w,j=3:2:h]);vec([(i,j) for i=3:2:w,j=2:2:h])]
        update!(cs, rule, idx[1], idx[2])
    end
end

function update_random!(cs::CellSpace, rule::Rule)
    for idx in shuffle(vec([(i,j) for i=2:cs.width-1,j=2:cs.height-1]))
        update!(cs, rule, idx[1], idx[2])
    end
end

function update!(cs::CellSpace, rule::Rule, x::Integer, y::Integer)
    # get 3x3 neighbor cells array
    ncells = view(cs.cells, x-1:x+1, y-1:y+1)

    # get current state of (x,y)
    state = get_target_state(ncells)

    # transition
    state,flag = transition(state, rule)

    # set state
    set_target_state(ncells, state)

    flag
end

function transition(key::Transition, rule)
    if haskey(rule.dict, key)
        rule.dict[key][1], true
    else
        key, false
    end
end

function save(fname::AbstractString, cs::CellSpace)
    cells = map(tostr, cs.cells)
    writedlm(fname, cells, ' ')
end

function load_cell(fname::AbstractString)
    tcells = readdlm(fname,String)
    load_cell(tcells)
end

function load_cell(tcells::Array{String,2})
    if !reduce(&, map(x->ismatch(r"^[0-9a-zA-Z]{4,4}$",x), tcells))
        throw("cell data error.")
    end
    cells = map(toval32, tcells)
    CellSpace(cells)
end

function get_state(cs::CellSpace, x::Integer, y::Integer, typ::Symbol)
    # typ is expected to be :north, :souch, :east, and :west
    eval(typ)(cs.cells[x, y])
end

function get_state(cs::CellSpace, x::Integer, y::Integer)
    cs.cells[x, y]
end

function show(io::IO, cs::CellSpace)
    cells = map(tostr, cs.cells)
    writedlm(io, cells, ' ')
end

function print(cs::CellSpace)
    show(STDOUT, cs)
end

function println(cs::CellSpace)
    show(STDOUT, cs)
    println()
end


text:
module KernelGradients

import Enzyme
import KernelAbstractions: Kernel

function Enzyme.autodiff(kernel::Kernel{<:Any, <:Any, <:Any, Fun}) where Fun
    f = kernel.f
    function df(ctx, args...)
        Enzyme.autodiff_deferred(f::Fun, Enzyme.Const, ctx, args...)
    end
    similar(kernel, df)
end

end # module


text:
module IPCC

using CairoMakie, GeoJSON, GeoMakie, Proj4

"""
	hexagons(df,clv,ttl,colors)

```
df=IPCC_hexagons()
include("examples/IPCC.jl")
clv, ttl, colors=example_hexagons(df)
f=IPCC.hexagons(df,clv,ttl,colors)
save("f.png", f)
```
"""
function hexagons(df,clv,ttl,colors)

	cl=["⚫⚫⚫","⚫⚫","⚫","⚪"]

	set_theme!(theme_light())
	g = Figure()

	a = Axis(g[1, 1])

	h=25

	txt="Type of observed change \nin hot extremes"
	t = text!(a, txt,
			position = Point2f(0,h*0.9),
			align = (:left, :baseline),
			justification = :left,
			color = :black,
			textsize = 12)

	poly!(hexa(1,h*0.8,0.5), color = :lightcoral, strokecolor = :black, strokewidth = 1)
	text!(a,"Increase (41)",position = Point2f(2,h*0.78),textsize = 10)
	poly!(hexa(1,h*0.7,0.5), color = :lightblue, strokecolor = :black, strokewidth = 1)
	text!(a,"Decrease (0)",position = Point2f(2,h*0.68),textsize = 10)
#	poly!(hexa(1,h*0.6,0.5), color = Pattern("/"), strokecolor = :black, strokewidth = 1)
	poly!(hexa(1,h*0.6,0.5), color = :yellow, strokecolor = :black, strokewidth = 1)
	text!(a,"Low agreement in the \n type of change (2)",position = Point2f(2,h*0.58),textsize = 10)
	poly!(hexa(1,h*0.5,0.5), color = :lightgray, strokecolor = :black, strokewidth = 1)
	text!(a,"Limited data and/or \n literature (2)",position = Point2f(2,h*0.48),textsize = 10)

	txt="Confidence in human \ncontribution to the \nobserved change"
	t = text!(a, txt,
			position = Point2f(0,h*0.3),
			align = (:left, :baseline),
			justification = :left,
			color = :black,
			textsize = 12)

	text!(a,cl[1]*" High",position = Point2f(1,h*0.24),textsize = 10)
	text!(a,cl[2]*"  Medium",position = Point2f(1,h*0.20),textsize = 10)
	text!(a,cl[3]*"    Low due to limited agreement",position = Point2f(1,h*0.16),textsize = 10)
	text!(a,cl[4]*"    Low due to limited evidence",position = Point2f(1,h*0.12),textsize = 10)

	hidespines!(a)
	hidedecorations!(a)
	xlims!(0,10)
	ylims!(0,h)

	a = Axis(g[1, 2],title=ttl)
	for i in 1:size(df,1)
		if iseven(df[i,:y])
			x0=cos(pi/6)*df[i,:x]
			y0=df[i,:y]*(1/2+sin(pi/6)/2)
			poly!(hexa(x0,y0,0.5), color = colors[i], strokecolor = :black, strokewidth = 1)
		else
			x0=cos(pi/6)*(df[i,:x]-0.5)
			y0=(0.5+sin(pi/6)/2)*df[i,:y]
			poly!(hexa(x0,y0,0.5), color = colors[i], strokecolor = :black, strokewidth = 1)
		end

		text!(
			df[i,:acronym],
			position = Point2f(x0,y0),
			align = (:center, :baseline),
			color = :black,
			textsize = 10,
		)

		ll=["⚪?"]
		for l in 1:4
			sum(clv[l].==df[i,:acronym])>0 ? ll[1]=cl[l] : nothing
		end
		ll[1]!=="⚪?" ? println("could not find confidenve level for $(df[i,:acronym])") : nothing

		text!(
			ll[1],
			position = Point2f(x0,y0-0.2),
			align = (:center, :baseline),
			color = :black,
			textsize = 10,
		)

	end

	text!(a,"Type of observed change since the 1950s",position = Point2f(0,-5),
		textsize = 16,align = (:center, :baseline), justification = :center)

	hidespines!(a)
	hidedecorations!(a)
	xlims!(-7,8)
	ylims!(-6.0,4.0)

	colsize!(g.layout, 1, Relative(1/5))

	g
end

#position of the hexagon points if center is index x0,y0 and radius is r
function hexa(x0,y0,r)
	hx=Point2f[]
	for i=1:6
		push!(hx,Point2f(x0+r*cos(-.5pi+(i-1)*pi/3),y0+r*sin(-.5pi+(i-1)*pi/3)))
	end
	hx
end

##

function fig1a(dat,dat1,dat2)
	set_theme!(theme_black())
	fig1a = Figure(resolution = (900, 600))

	##

	fig1a_ax1 = Axis(fig1a[1,1],xticksvisible=false,xticklabelsvisible=false,yticks=[0.2, 1.0])
	band!([0.0,1.0],fill(dat2.t5[1],2),fill(dat2.t95[1],2) ; color = (:greenyellow, 0.5))

	ylims!(-0.5,2.0)
	hidespines!(fig1a_ax1,:t,:r,:b)

	##

	fig1a_ax2 = Axis(fig1a[1,2], title="Change in global surface temperature", ylabel="degree C")

	l1=lines!(dat.year,dat.temp; color = :dodgerblue)
	band!(dat.year, dat.t5, dat.t95; color = (:dodgerblue, 0.5))

	l2=lines!(dat1.year,dat1.temp; color = :red,linewidth=2.0)

	band!([1850.0,2000.0],fill(-0.5,2),fill(2.0,2) ; color = (:greenyellow, 0.25))

	ylims!(-0.5,2.0)
	axislegend(fig1a_ax2,[l1, l2],["reconstructed","observed"],position=:lt)

	##

	colsize!(fig1a.layout, 2, Relative(0.95))

	fig1a
end

##

function fig1b(dat_1b)
	set_theme!(theme_black())
	fig1b = Figure(resolution = (900, 600))
	fig1b_ax2 = Axis(fig1b[1,1], title="Change in global surface temperature", ylabel="degree C")

	f1b_l1=lines!(dat_1b.year,dat_1b.HNm; color = :violet)
	band!(dat_1b.year,dat_1b.HN5,dat_1b.HN95; color = (:violet, 0.5))

	f1b_l2=lines!(dat_1b.year,dat_1b.Nm; color = :dodgerblue)
	band!(dat_1b.year,dat_1b.N5,dat_1b.N95; color = (:dodgerblue, 0.5))

	f1b_l3=lines!(dat_1b.year,dat_1b.obs; color = :yellow, linewidth=2.0)
	axislegend(fig1b_ax2,[f1b_l1, f1b_l2, f1b_l3],["Human+Nature","Nature","observed"],position=:lt)

	fig1b
end

##

function fig2(dat,dat_b,dat_c)

	set_theme!(theme_black())
	f = Figure(resolution = (1200, 800))

	f_ax1 = Axis(f[1,1],xticksvisible=false,xticklabelsvisible=false,
		title="Observed warming \n 2010-2019 relative \n to 1850-1900 \n \n")
	barplot!([dat[1,2]],color=:gray70)
	errorbars!([1.0], [dat[1,2]], [dat[1,2]-dat[2,2]], [dat[3,2]-dat[1,2]], color = :white, whiskerwidth = 20, linewidth=4.0)
	xlims!(-0.5,2.5); ylims!(-1.0,2.0)

	ttl1a="Aggregated contributions to \n 2010-2019 warming relative to \n 1850-1900, assessed from \n attribution studies."
	xti1a=names(dat_b)[2:end]; n1a=length(xti1a); ye=:yellow; wh=:white;
	xpo1a=1:n1a; xco1a=[:orange,ye,ye,wh,wh]
	bco1a=[:darkorchid1,:darkorchid1,:deepskyblue,wh,wh]
	f_ax2 = Axis(f[1,2],title=ttl1a,xticks = (xpo1a,xti1a),xticklabelrotation=.25pi,xticklabelcolor=xco1a)
	barplot!([dat_b[1,i] for i in 2:size(dat_b,2)],color=bco1a)
	ylims!(-1.0,2.0)

	xs=collect(2:size(dat_b,2)) .- 1.0
	ys=[dat_b[1,i] for i in 2:size(dat_b,2)]
	lowerrors=[dat_b[2,i]-dat_b[1,i] for i in 2:size(dat_b,2)]
	higherrors=[dat_b[2,i]-dat_b[1,i] for i in 2:size(dat_b,2)]
	errorbars!(xs, ys, lowerrors, higherrors, color = :white, whiskerwidth = 20, linewidth=4.0)

	ttl1a="Contributions to 2010-2019 \n warming relative to 1850-1900, \n assessed from radiative forcing studies \n "
	xti1a=dat_c[:,1]; n1a=length(xti1a);
	xpo1a=collect(1:n1a); xco1a=fill(wh,n1a)
	bco1a=fill(wh,n1a)
	bco1a[findall(dat_c[!,2].>0)].=:darkorchid1
	bco1a[findall(dat_c[!,2].<0)].=:deepskyblue
	f_ax3 = Axis(f[1,3],title=ttl1a,xticks=(xpo1a,xti1a),xticklabelcolor="white",xticklabelrotation=.25pi)
	barplot!([dat_c[i,2] for i in 1:n1a], color=bco1a)
	ylims!(-1.0,2.0)

	xs=collect(1:size(dat_c,1))
	ys=[dat_c[i,2] for i in 1:size(dat_c,1)]
	lowerrors=[dat_c[i,2]-dat_c[i,3] for i in 1:size(dat_c,1)]
	higherrors=[dat_c[i,2]-dat_c[i,3] for i in 1:size(dat_c,1)]
	errorbars!(xs, ys, lowerrors, higherrors, color = :white, whiskerwidth = 20, linewidth=4.0)

	colsize!(f.layout, 1, Relative(0.1))
	colsize!(f.layout, 2, Relative(0.2))
	colsize!(f.layout, 3, Relative(0.7))

	f
end

##

function fig4a(dat)
	nam=["SSP1-1.9","SSP1-2.6","SSP2-4.5","SSP3-7.0","SSP5-8.5"]
	var=["ssp119","ssp126","ssp245","ssp370","ssp585"]
	col=[:lightblue,:darkblue,:orange,:red,:darkred]

	g = Figure(resolution = (900, 600))
	a = Axis(g[1:3, 1],yticks=-20:20:140, title="Carbon dioxide (GtCO2/yr)")

	for i in 1:length(nam)
		lines!(dat.C.years,dat.C[!,var[i]],label=nam[i],color=col[i])
	end
	xlims!(2015,2100)

	Legend(g[2, 2], a)

	a = Axis(g[1, 3],yticks=0:200:800, title="Methane (MtCH4/yr)")
	for i in 1:length(nam)
		lines!(dat.M.years,dat.M[!,var[i]],label=nam[i],color=col[i])
	end
	xlims!(2015,2100); ylims!(0,800)

	a = Axis(g[2, 3],yticks=0:5:20, title="Nitrous oxide (MtN2O/yr)")
	for i in 1:length(nam)
		lines!(dat.N.years,dat.N[!,var[i]],label=nam[i],color=col[i])
	end
	xlims!(2015,2100); ylims!(0,20)

	a = Axis(g[3, 3],yticks=0:40:120, title="Sulfur dioxide (MtSO2/yr)")
	for i in 1:length(nam)
		lines!(dat.S.years,dat.S[!,var[i]],label=nam[i],color=col[i])
	end
	xlims!(2015,2100); ylims!(0,120)

	g
end

function mybarplot!(dat_b,nam="ssp126",col=:royalblue3,col2=:midnightblue)
	tmp1=dat_b[findall(dat_b.scenario.==nam),:]
	t=tmp1[findall(tmp1.forcing.=="total")[1],:]
	c=tmp1[findall(tmp1.forcing.=="co2")[1],:]
	n=tmp1[findall(tmp1.forcing.=="non_CO2_GHGs")[1],:]
	a=tmp1[findall(tmp1.forcing.=="Aersols_Landuse")[1],:]

	barplot!([t.p50,c.p50,n.p50,a.p50],color=col)

	o=dat_b[findall(dat_b.scenario.=="observations")[1],:]
	barplot!([o.p50],color=col2)

	errorbars!(collect(1.0:4.0),[t.p50,c.p50,n.p50,a.p50],
		[t.p50,c.p50,n.p50,a.p50].-[t.p05,c.p05,n.p05,a.p05],
		[t.p95,c.p95,n.p95,a.p95].-[t.p50,c.p50,n.p50,a.p50],
		color = :black, whiskerwidth = 20, linewidth=2.0)

	ylims!(-1.0,6.0)
end

function fig4b(dat_b)
	nam=["SSP1-1.9","SSP1-2.6","SSP2-4.5","SSP3-7.0","SSP5-8.5"]
	var=["ssp119","ssp126","ssp245","ssp370","ssp585"]

	set_theme!(theme_light())
	h = Figure(resolution = (900, 600))

	xtla=["total (observed)","CO2","Non-CO2 GHGs","Aerosols, land use"]

	b = Axis(h[1, 1],xticks=(1:4,xtla), title=nam[1], textsize=12, xticklabelrotation=.35pi)
	mybarplot!(dat_b,var[1],:deepskyblue,:deepskyblue3)

	b = Axis(h[1, 2],xticks=(1:4,xtla), title=nam[2], textsize=12, xticklabelrotation=.35pi)
	mybarplot!(dat_b,var[2],:royalblue3,:midnightblue)

	b = Axis(h[1, 3],xticks=(1:4,xtla), title=nam[3], textsize=12, xticklabelrotation=.35pi)
	mybarplot!(dat_b,var[3],:tan1,:tan3)

	b = Axis(h[1, 4],xticks=(1:4,xtla), title=nam[4], textsize=12, xticklabelrotation=.35pi)
	mybarplot!(dat_b,var[4],:firebrick2,:firebrick3)

	b = Axis(h[1, 5],xticks=(1:4,xtla), title=nam[5], textsize=12, xticklabelrotation=.35pi)
	mybarplot!(dat_b,var[5],:brown3,:brown4)

	h
end

##

function myproj(dat)
	source=Proj4.Projection("+proj=longlat +datum=WGS84")
	dest=Proj4.Projection("+proj=eqearth +lon_0=200.0 +lat_1=0.0 +x_0=0.0 +y_0=0.0 +ellps=GRS80")

	xy=Proj4.transform(source, dest, [vec(dat.lon) vec(dat.lat)])
	x=xy[:,1]
	y=xy[:,2]

	x=reshape(x,size(dat.lon))
	y=reshape(y,size(dat.lon))

	x=circshift(x, (-20,0))
	y=circshift(y, (-20,0))
	z=circshift(dat.var, (-20,0))

	(x=x,y=y,z=z)
end

function get_land_geo()
	url = "https://raw.githubusercontent.com/PublicaMundi/MappingAPI/master/data/geojson/"
	#https://raw.githubusercontent.com/nvkelso/natural-earth-vector/master/geojson/countries.geojson
	land = Downloads.download(url * "countries.geojson", IOBuffer())
    land_geo = GeoJSON.read(seekstart(land))
	land_geo.features[1]
end

#code adapted from https://lazarusa.github.io/BeautifulMakie/
function fig5(dat)
	projection_choice=1

	projection_choice==1 ? dx=-Int(size(dat.lon,1)/2) : dx=-20
	lons = circshift(dat.lon[:,1],dx)
	lats = dat.lat[1,:]
	field = circshift(dat.var,(dx,0))

	txt_source="+proj=longlat +datum=WGS84"
	txt_dest="+proj=eqearth +lon_0=200.0 +lat_1=0.0 +x_0=0.0 +y_0=0.0 +ellps=GRS80"
	if projection_choice==2
		trans = Proj4.Transformation(txt_source,txt_dest, always_xy=true)
	else
		trans = Proj4.Transformation("+proj=longlat +datum=WGS84", "+proj=wintri", always_xy=true)
	end
	#+proj=wintri, natearth2

	ptrans = Makie.PointTrans{2}(trans)
	fig = Figure(resolution = (1200,800), fontsize = 22)
	ax = Axis(fig[1,1], aspect = DataAspect(), title = dat.meta.ttl)
	# all input data coordinates are projected using this function
	ax.scene.transformation.transform_func[] = ptrans
	# add some limits, still it needs to be manual
	points = [Point2f0(lon, lat) for lon in lons, lat in lats]
	rectLimits = FRect2D(Makie.apply_transform(ptrans, points))
	limits!(ax, rectLimits)

	hm1 = surface!(ax, lons, lats, field, shading = false, overdraw = false,
	colorrange=dat.meta.colorrange, colormap=dat.meta.cmap)

	hm2 = lines!(ax, GeoMakie.coastlines(), color = :black, overdraw = true)

	##
	if projection_choice==1
		lonrange = -180:60:180
	else
		lonrange = collect(20.0:60:380); lonrange[end]-=0.0001
	end
	latrange = -90.0:30:90

	lonlines = [Point2f0(j,i) for i in lats, j in lonrange]
	latlines = [Point2f0(j,i) for j in lons, i in latrange]

	[lines!(ax, lonlines[:,i], color = (:black,0.25),
	 linestyle = :dash, overdraw = true) for i in 1:size(lonlines)[2]]
	[lines!(ax, latlines[:,i], color = (:black,0.25), linestyle = :dash,
	 overdraw = true) for i in 1:size(latlines)[2]]

	xticks = first.(trans.(Point2f0.(lonrange, -90)))
	yticks = last.(trans.(Point2f0.(-180,latrange)))
	ax.xticks = (xticks, string.(lonrange, 'ᵒ'))
	ax.yticks = (yticks, string.(latrange, 'ᵒ'))

	#add colorbar
	Colorbar(fig[1,2], hm1, height = Relative(0.65))

	# hide just original grid
	hidedecorations!(ax, ticks = false, label = false, ticklabels=false)
	hidespines!(ax)

	##

	#save(joinpath(@__DIR__, "change_at_1C.png"), fig) # HIDE

	fig
end

end #module IPCC


text:
using HJBSolver, FactCheck

let

    function createmodel()
        T = 3.
        γ = 5e-2
        C = -1.

        xmin = 0.; xmax = 1.
        αmin = 0.; αmax = 1.

        d(α) = 1.-α
        b(t,x,α) = -d(α)
        σ(t,x,α) = γ*d(α)
        f(t,x,α) = α*d(α)
        g(x) = C*x
        Dmin(t) = zero(t)

        model = HJBOneDim(b, σ, f, g, T, αmin, αmax, xmin, xmax,
                          (true, false), (Dmin, zero))
    end

    facts("1D, no boundary condition, Constant Policy") do
        K = 101; N = 60;  M = 80
        model = createmodel()

        @time v, pol = solve(model, K, N, M)

        vcheck = 0.6230938105907916
        polcheck = 0.6455696202531646

        @fact v[end,1] --> roughly(vcheck, 1e-5)
        @fact pol[end,1] --> roughly(polcheck, 1e-5)
    end

    facts("1D, no boundary condition, Policy-iteration") do
        K = 101; N = 60;
        model = createmodel()

        @time v, pol = solve(model, K, N)

        vcheck = 0.6246558092538466
        polcheck = 0.6398206678262762

        @fact v[end,1] --> roughly(vcheck, 1e-5)
        @fact pol[end,1] --> roughly(polcheck, 1e-5)
    end

end

# TODO: do it for opposite boundary as well

max_stars_repo_name: anriseth/HJBSolver

max_issues_repo_name: anriseth/HJBSolver

max_forks_repo_name: anriseth/HJBSolver


text:
"
Auxiliary structure that keeps some generic information for GEP experiments.

# Fields
- `_functionSet::Array{FunctionNode}`: array that contains all the non-terminal
    nodes specified in the node information file.
- `_terminalSet::Array{TerminalNode}`: array that contains all the terminal nodes
    specified in the node information file.
- `_head::UInt16`: length of the head of the genotype.
- `_tail::UInt16`: length of the tail of the genotype. It is calculated as follows:
    t = h (n-1) + 1 , where n is the maximum arity of the non-terminal nodes of the
    problem.
- `_variables::Array{String}`: array that contains the names of the variables of
    the problem.

See also: [`GPExperimentInfo`](@ref)
"
struct GEPInfo <: GPExperimentInfo
    _functionSet::Array{FunctionNode}
    _terminalSet::Array{TerminalNode}
    _head::UInt16
    _tail::UInt16
    _variables::Array{String}
end # struct

max_stars_repo_name: SergioGmezM/GeneticJulia

max_issues_repo_name: SergioGmezM/GeneticJulia

max_forks_repo_name: SergioGmezM/GeneticJulia


text:
@testset "RockPaperScissorsEnv" begin

    rng = StableRNG(123)
    env = RockPaperScissorsEnv()

    RLBase.test_interfaces!(env)
    RLBase.test_runnable!(env)

    rewards = [[], []]
    for _ in 1:50_000
        while !is_terminated(env)
            env(rand(rng, legal_action_space(env)))
        end
        push!(rewards[1], reward(env, 1))
        push!(rewards[2], reward(env, 2))
        reset!(env)
    end

    @test all(rewards[1] .+ rewards[2] .== 0)

end


text:
using SpecialFunctions
import Base.Broadcast

const constant_funcs = []

const monadic_linear = [deg2rad, +, rad2deg, transpose, -, conj]

const monadic_nonlinear = [asind, log1p, acsch, erfc, digamma, acos, asec, acosh, airybiprime, acsc, cscd, log, tand, log10, csch, asinh, airyai, abs2, gamma, lgamma, erfcx, bessely0, cosh, sin, cos, atan, cospi, cbrt, acosd, bessely1, acoth, erfcinv, erf, dawson, inv, acotd, airyaiprime, erfinv, trigamma, asecd, besselj1, exp, acot, sqrt, sind, sinpi, asech, log2, tan, invdigamma, airybi, exp10, sech, erfi, coth, asin, cotd, cosd, sinh, abs, besselj0, csc, tanh, secd, atand, sec, acscd, cot, exp2, expm1, atanh]

diadic_of_linearity(::Val{(true, true, true)}) = [+, rem2pi, -, >, isless, <, isequal, max, min, convert]
diadic_of_linearity(::Val{(true, true, false)}) = [*]
#diadic_of_linearit(::(Val{(true, false, true)}) = [besselk, hankelh2, bessely, besselj, besseli, polygamma, hankelh1]
diadic_of_linearity(::Val{(true, false, false)}) = [/]
diadic_of_linearity(::Val{(false, true, false)}) = [\]
diadic_of_linearity(::Val{(false, false, false)}) = [hypot, atan, mod, rem, lbeta, ^, beta]
diadic_of_linearity(::Val) = []

const all_funcs = vcat(monadic_linear,
                       monadic_nonlinear,
                       diadic_of_linearity(Val{(true, false, false)}()),
                       diadic_of_linearity(Val{(false, true, false)}()),
                       diadic_of_linearity(Val{(false, false, false)}())) |> unique

isleaf(::Any) = false
isleaf(t::Tagged) = isleaf(t.value)
for f in all_funcs
    @eval isleaf(::$(typeof(f))) = true
end

haslinearity(f, nargs) = false

# some functions strip the linearity metadata

for f in constant_funcs
    @eval begin
        haslinearity(::typeof($f), ::Val) = true
        linearity(::typeof($f), ::Val) = nothing
    end
end

# linearity of a single input function is either
# Val{true}() or Val{false}()
#
for f in monadic_linear
    @eval begin
        haslinearity(::typeof($f), ::Val{1}) = true
        linearity(::typeof($f), ::Val{1}) = Val{true}()
    end
end
# linearity of a 2-arg function is:
# Val{(linear11, linear22, linear12)}()
#
# linearIJ refers to the zeroness of d^2/dxIxJ
for f in monadic_nonlinear
    @eval begin
        haslinearity(::typeof($f), ::Val{1}) = true
        linearity(::typeof($f), ::Val{1}) = Val{false}()
    end
end

for linearity_mask = 0:2^3-1
    lin = Val{map(x->x!=0, (linearity_mask & 4,
                            linearity_mask & 2,
                            linearity_mask & 1))}()

    for f in diadic_of_linearity(lin)
        @eval begin
            haslinearity(::typeof($f), ::Val{2}) = true
            linearity(::typeof($f), ::Val{2}) = $lin
        end
    end
end


text:
using Tables:matrix

"""
    UnsupervisedDetector

This abstract type forms the basis for all implemented unsupervised outlier detection algorithms. To implement a new
`UnsupervisedDetector` yourself, you have to implement the `fit(detector, X)::Fit` and
`score(detector, model, X)::Result` methods.
"""
abstract type UnsupervisedDetector <: MMI.Unsupervised end

"""
    SupervisedDetector

This abstract type forms the basis for all implemented supervised outlier detection algorithms. To implement a new
`SupervisedDetector` yourself, you have to implement the `fit(detector, X, y)::Fit` and
`score(detector, model, X)::Result` methods.
"""
abstract type SupervisedDetector <: MMI.Deterministic end

"""
    Detector::Union{<:SupervisedDetector, <:UnsupervisedDetector}

The union type of all implemented detectors, including supervised, semi-supervised and unsupervised detectors. *Note:* A
semi-supervised detector can be seen as a supervised detector with a specific class representing unlabeled data.
"""
const Detector = Union{<:SupervisedDetector,<:UnsupervisedDetector}

"""
    Evaluator

A evaluator uses one or more detector results produced by [`score`](@ref) and transforms them into labels or scores,
typically with two classes describing inliers (`1`) and outliers (`-1`).
"""
abstract type Evaluator <: MMI.Static end

"""
    Model

A `Model` represents the learned behaviour for specific [`Detector`](@ref). This might include parameters in parametric
models or other repesentations of the learned data in nonparametric models. In essence, it includes everything required
to transform an instance to an outlier score.
"""
abstract type Model end

"""
    Scores::AbstractVector{<:Real}

Scores are continuous values, where the range depends on the specific detector yielding the scores. *Note:* All
detectors return increasing scores and higher scores are associated with higher outlierness. Concretely, scores are
defined as an `AbstractVector{<:Real}`.
"""
const Scores = AbstractVector{<:Real}

"""
    Data::AbstractArray{<:Real}

The raw input data for every detector is defined as`AbstractArray{<:Real}` and should be a one observation per column
n-dimensional array. The input data used to [`fit`](@ref) a [`Detector`](@ref) and [`score`](@ref) [`Data`](@ref).
"""
const Data = AbstractArray{<:Real}

"""
    Labels::AbstractVector{<:Integer}

Labels are used for supervision and evaluation and are defined as an `AbstractArray{<:Integer}`. The convention for
labels is that `-1` indicates outliers, `1` indicates inliers and `0` indicates unlabeled data in semi-supervised tasks.
"""
const Labels = AbstractVector{<:Integer}

"""
    Fit

A `Fit` bundles a [`Model`](@ref) and [`Scores`](@ref) achieved when fitting a [`Detector`](@ref). The model
is used directly in later [`score`](@ref) calls and the (train-) scores are forwarded in [`score`](@ref).
"""
struct Fit
    model::Model
    scores::Scores
end

"""
    Result::Tuple{Scores, Scores}

Describes the result of using [`score`](@ref) with a [`Detector`](@ref) and is a tuple containing the train scores
achieved with [`fit`](@ref) and the test scores achieved with [`score`](@ref). We return both train and test scores
because this gives us the greatest flexibility in later score combination or classification.
"""
const Result = Tuple{Scores,Scores}

const _input_data = """    X::Union{AbstractMatrix, Tables.jl-compatible}
Either a matrix or a [Tables.jl-compatible](https://github.com/JuliaData/Tables.jl) data source, with one observation
per row and a number of feature columns."""

const _label_data = """    y::AbstractVector{<:Integer}
A vector of labels with `-1` indicating an outlier and `1` indicating an inlier."""

const _detector = """    detector::Detector
Any [`UnsupervisedDetector`](@ref) or [`SupervisedDetector`](@ref) implementation."""

_score_unsupervised(name::String) = """
```julia
using OutlierDetection: $name, fit, score
detector = $name()
X = rand(10, 100)
model = fit(detector, X)
train_scores, test_scores = score(detector, model, X)
```"""

_score_supervised(name::String) = """
```julia
using OutlierDetection: $name, fit, score
detector = $name()
X = rand(10, 100)
y = rand([-1,1], 100)
model = fit(detector, X, y)
train_scores, test_scores = score(detector, model, X)
```"""

const _evaluator = """
```julia
using OutlierDetection: Class, KNN, fit, score
detector = KNN()
X = rand(10, 100)
model = fit(detector, X)
train_scores, test_scores = score(detector, model, X)
ŷ = detect(Class(), train_scores, test_scores)
# or, if using multiple detectors
# ŷ = detect(Class(), (train1, test1), (train2, test2), ...)
```"""

"""
    fit(detector,
        X,
        y)

Fit a specified unsupervised, supervised or semi-supervised outlier detector. That is, learn a `Model` from input data
`X` and, in the supervised and semi-supervised setting, labels `y`. In a supervised setting, the label `-1` represents
outliers and `1` inliers. In a semi-supervised setting, the label `0` additionally represents unlabeled data. *Note:*
Unsupervised detectors can be fitted without specifying `y`, otherwise `y` is simply ignore.

Parameters
----------
$_detector

$_input_data

Returns
----------
    fit::Fit
The learned model of the given detector, which contains all the necessary information for later prediction and the
achieved outlier scores of the given input data `X`.

Examples
--------
$(_score_unsupervised("KNN"))
""" # those definitions apply when X is not already a (transposed) abstract array
fit(detector::UnsupervisedDetector, X) = fit(detector, matrix(X; transpose=true)) # unsupervised call syntax
fit(detector::SupervisedDetector, X, y) = fit(detector, matrix(X; transpose=true), y)

"""
    score(detector,
          model,
          X)

Transform input data `X` to outlier scores using an [`UnsupervisedDetector`](@ref) or [`SupervisedDetector`](@ref) and
a corresponding [`Model`](@ref).

Parameters
----------
$_detector

    model::Model
The model learned from using [`fit`](@ref) with a supervised or unsupervised [`Detector`](@ref)

$_input_data

Returns
----------
    result::Result
Tuple of the achieved outlier scores of the given train and test data.

Examples
--------
$(_score_unsupervised("KNN"))
""" # definition applies when X is not already a (transposed) abstract array
score(detector::Detector, fitresult::Fit, X) = score(detector, fitresult, matrix(X; transpose=true))

"""
    detect(Evaluator,
           result...)

Convert a number of scores into inlier (`1`) / outlier (`-1`) classes, typically by using on a outlier-threshold on the
achieved scores.

Parameters
----------
    evaluator::Evaluator
A [`Evaluator`](@ref) that implements the [`detect`](@ref) method.

    result::Result...
One or more [`score`](@ref) results (tuples) or alternatively a single vector of scores or two vectors, where the first
vector represents train scores and the second vector test scores.

Returns
----------
    result::Labels
A vector containing the binary inlier and outlier labels.

Examples
--------
$_evaluator
""" # transforms single scores, or combination of train-test scores into a tuple
detect(evaluator::Evaluator, scores_train::Scores) = detect(evaluator, (scores_train, scores_train))
detect(evaluator::Evaluator, scores_train::Scores, scores_test::Scores) =
    detect(evaluator, (scores_train, scores_test))

"""
    @score

Helps with the definition of [`score`](@ref) for detectors, by unpacking the `.model` field of the second argument
directly into the argument and implicitly returning the values of the `.scores` field as the first tuple element of the
returned expression.
"""
macro score(fn)
    fn = MacroTools.longdef(fn)
    @capture(fn, function f_(detector_, result_::Fit, X_::Data)::Result body_ end) || error("Expected a function with
    three parameters f(`detector<:Detector`, `result::Fit`, `X::Data`)::Result and fully specified types.")
    copy_result = :copy_result_unique_name

    # implicitly return the scores as the first tuple element for all return
    return_count = 0
    body = MacroTools.postwalk(body) do x
        @capture(x, ret_return) || return x
        # count the return
        return_count += 1

        # get the return expression
        ret_expr = ret.args[1]
        return :(($copy_result.scores, $(ret_expr)))
    end

    # check if there are any returns were found, if not, use the last body expression
    if return_count == 0
        body = :(($copy_result.scores, $(body.args[end])))
    end

    :(function $f($detector, $result::Fit, $X::Data)::Result
        $copy_result = $result;
        $result = $result.model;
        $body
    end)
end


text:
# module for sensing object and drawing observation

using Plots, Random, Distributions
pyplot()

include(joinpath(split(@__FILE__, "src")[1], "src/model/map/map.jl"))
include(joinpath(split(@__FILE__, "src")[1], "src/common/observation_function/observation_function.jl"))

mutable struct Sensor
  map
  last_data
  dist_rng
  dir_rng
  dist_noise_rate
  dir_noise # std dev
  dist_bias_rate_std
  dir_bias
  phantom_dist_x
  phantom_dist_y
  phantom_prob
  oversight_prob
  occlusion_prob

  # init
  function Sensor(map::Map;
                  dist_rng::Tuple=(0.5, 6.0),
                  dir_rng::Tuple=(-pi/3, pi/3),
                  dist_noise_rate=0.0,
                  dir_noise=0.0,
                  dist_bias_rate_stddev=0.0,
                  dir_bias_stddev=0.0,
                  phantom_prob=0.0,
                  phantom_rng_x::Tuple=(-5.0, 5.0),
                  phantom_rng_y::Tuple=(-5.0, 5.0),
                  oversight_prob=0.0,
                  occlusion_prob=0.0)
    self = new()
    self.map = map
    self.last_data = []
    self.dist_rng = dist_rng
    self.dir_rng = dir_rng
    self.dist_noise_rate = dist_noise_rate
    self.dir_noise = dir_noise
    self.dist_bias_rate_std = rand(Normal(0.0, dist_bias_rate_stddev))
    self.dir_bias = rand(Normal(0.0, dir_bias_stddev))
    rx = phantom_rng_x
    ry = phantom_rng_y
    self.phantom_dist_x = Uniform(rx[1], rx[2])
    self.phantom_dist_y = Uniform(ry[1], ry[2])
    self.phantom_prob = phantom_prob
    self.oversight_prob = oversight_prob
    self.occlusion_prob = occlusion_prob
    return self
  end
end

function visible(self::Sensor, obsrv=nothing)
  if obsrv === nothing
    return false
  end

  return ((self.dist_rng[1] <= obsrv[1] <= self.dist_rng[2])
    && (self.dir_rng[1] <= obsrv[2] <= self.dir_rng[2]))
end

function noise(self::Sensor, obsrv::Array)
  ell = rand(Normal(obsrv[1], obsrv[1] * self.dist_noise_rate))
  phi = rand(Normal(obsrv[2], self.dir_noise))
  return [ell, phi]
end

function bias(self::Sensor, obsrv::Array)
  return obsrv + [obsrv[1] * self.dist_bias_rate_std, self.dir_bias]
end

function phantom(self::Sensor, sns_pose::Array, obsrv::Array)
  if rand(Uniform()) < self.phantom_prob
    phantom_pose = [rand(self.phantom_dist_x), rand(self.phantom_dist_y)]
    return observation_function(sns_pose, phantom_pose)
  else
      return obsrv
  end
end

function oversight(self::Sensor, obsrv::Array)
  if rand(Uniform()) < self.oversight_prob
    return nothing
  else
    return obsrv
  end
end

function occlusion(self::Sensor, obsrv::Array)
  if rand(Uniform()) < self.occlusion_prob
    ell = obsrv[1] - rand(Uniform()) * (self.dist_rng[2] - obsrv[1])
    phi = obsrv[2]
    return [ell, phi]
  else
    return obsrv
  end
end

function data(self::Sensor, sns_pose::Array)
  observed = []
  for obj in self.map.objects
    obsrv = observation_function(sns_pose, obj.pose)
    obsrv = phantom(self, sns_pose, obsrv)
    obsrv = occlusion(self, obsrv)
    obsrv = oversight(self, obsrv)
    if visible(self, obsrv)
      obsrv = bias(self, obsrv)
      obsrv = noise(self, obsrv)
      push!(observed, (obsrv, obj.id))
    end
  end
  self.last_data = observed
end

function draw!(self::Sensor, sns_pose::Array)
  for obj in self.last_data
    x, y, theta = sns_pose[1], sns_pose[2], sns_pose[3]
    distance, direction = obj[1][1], obj[1][2]
    ox = x + distance * cos(direction + theta)
    oy = y + distance * sin(direction + theta)
    plot!([x, ox], [y, oy], color="pink",
          legend=false, aspect_ratio=true)
  end
end

max_stars_repo_name: ShisatoYano/JuliaAutonomy

max_issues_repo_name: ShisatoYano/JuliaAutonomy

max_forks_repo_name: ShisatoYano/JuliaAutonomy


text:
# Many derivative free methods of different orders
#
# TODO: rework Order5 #https://pdfs.semanticscholar.org/ce50/3210d96f653a14b28da96600d5990d2abe97.pdf
# https://content.sciendo.com/view/journals/tmj/10/4/article-p103.xml 7 and 8
# order8: https://www.hindawi.com/journals/ijmms/2012/493456/ref/

## Order0 and 1 are secant type
function init_state(method::AbstractSecant, fs, x::Number)
    x1 = float(x)
    x0 = _default_secant_step(x1)
    init_state(method, fs, (x0, x1))
end

function init_state(method::AbstractSecant, fs, x::Union{Tuple, Vector})
    x0, x1 = promote(float(x[1]), float(x[2]))
    fx0, fx1 = fs(x0), fs(x1)
    state = UnivariateZeroState(x1, x0, eltype(x1)[],
                                fx1, fx0, eltype(fx1)[],
                                0, 2,
                                false, false, false, false, "")

    state
end

function init_state!(state::UnivariateZeroState{T, S}, method::AbstractSecant, fs, x::Number) where {T, S}
    x1::T = float(x)
    x0::T = _default_secant_step(x1)
    init_state!(state, method, fs, (x0, x1))
end

function init_state!(state::UnivariateZeroState{T, S}, ::AbstractSecant, f, x::Union{Tuple, Vector}) where {T, S}
    x0,x1 = promote(float.(x)...)
    fx0, fx1 = promote(f(x0), f(x1))
    init_state!(state, x1, x0, T[], fx1, fx0, S[])
    state.fnevals = 2
    nothing
end

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

## Order0 and Secant are related
"""
    Order0()

The `Order0` method is engineered to be a more robust, though possibly
slower, alternative to to the other derivative-free root-finding
methods. The implementation roughly follows the algorithm described in
*Personal Calculator Has Key to Solve Any Equation f(x) = 0*, the
SOLVE button from the
[HP-34C](http://www.hpl.hp.com/hpjournal/pdfs/IssuePDFs/1979-12.pdf).
The basic idea is to use a secant step. If along the way a bracket is
found, switch to bisection, using `AlefeldPotraShi`.  If the secant
step fails to decrease the function value, a quadratic step is used up
to 3 times.

"""
struct Order0 <: AbstractSecant end

function find_zero(fs, x0, method::Order0;
                   tracks::AbstractTracks=NullTracks(),
                   verbose=false,
                   kwargs...)
    M = Order1()
    N = AlefeldPotraShi()
    _find_zero(fs, x0, M, N; tracks=tracks,verbose=verbose, kwargs...)
end

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

## Secant
## https://en.wikipedia.org/wiki/Secant_method
"""
    Order1()

The `Order1()` method is an alias for `Secant`. It specifies the
[secant method](https://en.wikipedia.org/wiki/Secant_method).
This method keeps two values in its state, `x_n` and `x_n1`. The
updated point is the intersection point of x axis with the secant line
formed from the two points. The secant method uses 1 function
evaluation per step and has order `(1+sqrt(5))/2`.

"""
struct Secant <: AbstractSecant end
const Order1 = Secant

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

    if (o.fxn0 == o.fxn1) || (o.xn0 == o.xn1)
         o.stopped = true
         o.message = "Derivative approximation had issues"
         return
     end

    dx = o.fxn1 * (o.xn1 - o.xn0) / (o.fxn1 - o.fxn0)
    o.xn0, o.xn1 = o.xn1, o.xn1 - dx
    o.fxn0, o.fxn1 = o.fxn1, fs(o.xn1)
    incfn(o)

    nothing

end

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

### Steffensen
## https://en.wikipedia.org/wiki/Steffensen's_method#Simple_description
struct Steffensen <: AbstractUnivariateZeroMethod
end

"""
    Order2()

The quadratically converging
[Steffensen](https://en.wikipedia.org/wiki/Steffensen's_method#Simple_description)
method is used for the derivative-free `Order2()` algorithm. Unlike
the quadratically converging Newton's method, no derivative is
necessary, though like Newton's method, two function calls per step
are. This algorithm is more sensitive than Newton's method to poor
initial guesses.

"""
const Order2 = Steffensen

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

    wn::T = o.xn1 + steff_step(o.xn1, o.fxn1)

    fwn::S = fs(wn)
    incfn(o)

    fp, issue = _fbracket(o.xn1, wn, o.fxn1, fwn)

    if issue
        o.stopped = true
        o.message = "Derivative approximation had issues"
        return
    end

    o.xn0 = o.xn1
    o.fxn0 = o.fxn1
    o.xn1 = o.xn1 - o.fxn1 / fp #xn1
    o.fxn1 = fs(o.xn1)
    incfn(o)

    nothing
end

steffenson(f, x0; kwargs...) = find_zero(f, x0, Steffensen(); kwargs...)

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

"""
    Order5()

Implements an order 5 algorithm from *A New Fifth Order Derivative
Free Newton-Type Method for Solving Nonlinear Equations* by Manoj
Kumar, Akhilesh Kumar Singh, and Akanksha, Appl. Math. Inf. Sci. 9,
No. 3, 1507-1513 (2015), DOI: 10.12785/amis/090346. Four function
calls per step are needed.

"""
struct Order5 <: AbstractUnivariateZeroMethod end

## If we have a derivative, we have this
function update_state(method::Order5, fs::Union{FirstDerivative,SecondDerivative},
                      o::UnivariateZeroState{T,S}, options)  where {T, S}

    xn, fxn = o.xn1, o.fxn1

    fpxn::S = fs(xn, 1)
    incfn(o)

    if isissue(fpxn)
        o.stopped  = true
        return
    end

    yn::T = xn - fxn / fpxn
    fyn::S, fpyn::S = fs(yn), fs(yn, 1)
    incfn(o, 2)

    if isissue(fpyn)
        o.xn0, o.xn1 = xn, yn
        o.fxn0, o.fxn1 = fxn, fyn
        o.stopped  = true
        return
    end

    zn::T = xn  - (fxn + fyn) / fpxn
    fzn::S = fs(zn)
    incfn(o, 1)

    xn1 = zn - fzn / fpyn
    fxn1 = fs(xn1)
    incfn(o, 1)

    o.xn0, o.xn1 = xn, xn1
    o.fxn0, o.fxn1 = fxn, fxn1

    nothing
end

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

    xn = o.xn1
    fxn = o.fxn1

    wn::T = o.xn1 + steff_step(o.xn1, o.fxn1)

    fwn::S = fs(wn)
    incfn(o)

    fp, issue = _fbracket(o.xn1, wn, o.fxn1, fwn)
    if issue
        o.xn0, o.xn1 = o.xn1, wn
        o.fxn0, o.fxn1 = o.fxn1, fwn
        o.message = "Issue with divided difference f[xn, wn]"
        o.stopped  = true
        return
    end

    yn::T = o.xn1 - o.fxn1 / fp
    fyn::S = fs(yn)
    incfn(o)

    zn::T = xn - (fxn + fyn) / fp
    fzn::S = fs(zn)
    incfn(o)

    fp, issue = _fbracket_ratio(yn, o.xn1, wn, fyn, o.fxn1, fwn)
    if issue
        o.xn0, o.xn1 = o.xn1, yn
        o.fxn0, o.fxn1 = o.fxn1, fyn
        o.message = "Issue with f[xn,yn]*f[yn,wn] / f[xn, wn]"
        o.stopped = true
        return
    end

    o.xn0 = o.xn1
    o.fxn0 = o.fxn1
    o.xn1 = zn  - fzn  / fp
    o.fxn1 = fs(o.xn1)
    incfn(o)

    nothing
end

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

## cf also: https://doi.org/10.1515/tmj-2017-0049
"""
    Order8()

Implements an eighth-order algorithm from *New Eighth-Order
Derivative-Free Methods for Solving Nonlinear Equations* by Rajinder
Thukral, International Journal of Mathematics and Mathematical
Sciences Volume 2012 (2012), Article ID 493456, 12 pages DOI:
10.1155/2012/493456. Four function calls per step are required.

"""
struct Order8 <: AbstractUnivariateZeroMethod
end

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

    xn = o.xn1
    fxn = o.fxn1

    wn::T = xn + steff_step(xn, fxn)
    fwn::S = fs(wn)
    incfn(o)

    if isissue(fwn)
        o.xn0,o.xn1 = xn, wn
        o.fxn0,o.fxn1 = fxn, fwn
        o.stopped = true
        o.message = "issue with Steffensen step fwn"
        return
    end

    fp, issue = _fbracket(xn, wn, fxn, fwn)

    if issue
        o.stopped = true
        o.message = "issue with divided difference f[xn, wn]"
        return
    end

    yn::T = xn - fxn / fp
    fyn::S = fs(yn)
    incfn(o)

    fp, issue = _fbracket(yn, xn, fyn, fxn)
    if issue #fp
        o.xn0,o.xn1 = xn, yn
        o.fxn0,o.fxn1 = fxn, fyn
        o.stopped = true
        o.message = "issue with divided difference f[xn, yn]"
        return
    end

    phi = (1 + fyn / fwn)           # pick one of options
    zn =  yn - phi * fyn / fp
    fzn::S = fs(zn)
    incfn(o)

    fp, issue =  _fbracket_diff(xn, yn, zn, fxn, fyn, fzn)
    if issue
        o.xn0,o.xn1 = xn, zn
        o.fxn0,o.fxn1 = fxn, fzn
        o.message = "issue with divided difference  f[y,z] - f[x,y] + f[x,z]"
        o.stopped = true
        return
    end

    w = 1 / (1 - fzn/fwn)

    xi = (1 - 2fyn*fyn*fyn / (fwn * fwn * fxn))

    xn1 = zn - w * xi * fzn / fp
    fxn1::S = fs(xn1)
    incfn(o)

    o.xn0,o.xn1 = xn, xn1
    o.fxn0,o.fxn1 = fxn, fxn1

    nothing
end

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

"""
    Order16()

Implement the algorithm from
*New Sixteenth-Order Derivative-Free Methods for Solving Nonlinear Equations*
by R. Thukral,
American Journal of Computational and Applied Mathematics
p-ISSN: 2165-8935;    e-ISSN: 2165-8943; 2012;  2(3): 112-118
doi: 10.5923/j.ajcam.20120203.08.

Five function calls per step are required. Though rapidly converging,
this method generally isn't faster (fewer function calls/steps) over
other methods when using `Float64` values, but may be useful for
solving over `BigFloat`.

"""
struct Order16 <: AbstractUnivariateZeroMethod
end

function update_state(method::Order16, fs, o::UnivariateZeroState{T,S}, options) where {T, S}
    xn = o.xn1
    fxn = o.fxn1

    wn::T = xn + steff_step(xn, fxn)
    fwn::S = fs(wn)
    incfn(o)

    fp, issue = _fbracket(xn, wn, fxn, fwn)

    if issue
        o.xn0, o.xn1 = xn, wn
        o.fxn0, o.fxn1 = fxn, fwn
        o.message = "issue with f[xn,wn]"
        o.stopped = true
        return
    end

    yn::T = xn - fxn / fp
    fyn::S = fs(yn)
    incfn(o)

    fp, issue = _fbracket_ratio(yn, xn, wn, fyn, fxn, fwn)
    if issue
        o.xn0, o.xn1 = xn, yn
        o.fxn0, o.fxn1 = fxn, fyn
        o.message = "issue with f[xn,yn]*f[yn,wn]/f[xn,wn]"
        o.stopped = true
        return
    end

    zn = yn - fyn / fp
    fzn::S = fs(zn)
    incfn(o)

    fp, issue = _fbracket_diff(xn, yn, zn, fxn, fyn, fzn)
    u2, u3, u4 = fzn/fwn, fyn/fxn, fyn/fwn

    eta = 1 / (1 + 2*u3*u4^2) / (1 - u2)
    if issue
        o.xn0, o.xn1 = xn, zn
        o.fxn0, o.fxn1 = fxn, fzn
        o.stopped = true
        o.message = "Approximate derivative failed"
        return
    end

    an = zn - eta * fzn / fp
    fan::S = fs(an)
    incfn(o)

    fp, issue = _fbracket_ratio(an, yn, zn, fan, fyn, fzn)
    if issue
        o.xn0, o.xn1 = xn, an
        o.fxn0, o.fxn1 = fxn, fan
        o.stopped = true
        o.message = "Approximate derivative failed"
        return
    end

    u1, u5, u6 = fzn/fxn, fan/fxn, fan/fwn

    fp1, issue = _fbracket(xn,yn, fxn, fyn)

    sigma =  1 + u1*u2 - u1*u3*u4^2 + u5 + u6 + u1^2*u4 +
    u2^2*u3 + 3*u1*u4^2*(u3^2 - u4^2)/(fp1/oneunit(fp1))

    xn1 = an - sigma * fan / fp
    fxn1::S = fs(xn1)
    incfn(o)

    o.xn0, o.xn1 = xn, xn1
    o.fxn0, o.fxn1 = fxn, fxn1

    nothing
end


text:
type CurrentPlot
  nullableplot::Nullable{PlottingObject}
end
const CURRENT_PLOT = CurrentPlot(Nullable{PlottingObject}())

isplotnull() = isnull(CURRENT_PLOT.nullableplot)

function current()
  if isplotnull()
    error("No current plot/subplot")
  end
  get(CURRENT_PLOT.nullableplot)
end
current(plot::PlottingObject) = (CURRENT_PLOT.nullableplot = Nullable(plot))

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

Base.string(plt::Plot) = "Plot{$(plt.backend) n=$(plt.n)}"
Base.print(io::IO, plt::Plot) = print(io, string(plt))
Base.show(io::IO, plt::Plot) = print(io, string(plt))

getplot(plt::Plot) = plt
getplotargs(plt::Plot, idx::Int = 1) = plt.plotargs
convertSeriesIndex(plt::Plot, n::Int) = n

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

"""
The main plot command.  Use `plot` to create a new plot object, and `plot!` to add to an existing one:

```
  plot(args...; kw...)                  # creates a new plot window, and sets it to be the current
  plot!(args...; kw...)                 # adds to the `current`
  plot!(plotobj, args...; kw...)        # adds to the plot `plotobj`
```

There are lots of ways to pass in data, and lots of keyword arguments... just try it and it will likely work as expected.
When you pass in matrices, it splits by columns.  See the documentation for more info.
"""

# this creates a new plot with args/kw and sets it to be the current plot
function plot(args...; kw...)
  pkg = backend()
  d = Dict(kw)
  preprocessArgs!(d)
  dumpdict(d, "After plot preprocessing")

  plotargs = merge(d, getPlotArgs(pkg, d, 1))
  dumpdict(plotargs, "Plot args")
  plt = _create_plot(pkg; plotargs...)  # create a new, blank plot

  delete!(d, :background_color)
  plot!(plt, args...; d...)  # add to it
end

# this adds to the current plot, or creates a new plot if none are current
function  plot!(args...; kw...)
  local plt
  try
    plt = current()
  catch
    return plot(args...; kw...)
  end
  plot!(current(), args...; kw...)
end

# not allowed:
function plot!(subplt::Subplot, args...; kw...)
  error("Can't call plot! on a Subplot!")
end

# this adds to a specific plot... most plot commands will flow through here
function plot!(plt::Plot, args...; kw...)
  d = Dict(kw)
  preprocessArgs!(d)

  # for plotting recipes, swap out the args and update the parameter dictionary
  args = _apply_recipe(d, args...; kw...)

  dumpdict(d, "After plot! preprocessing")

  warnOnUnsupportedArgs(plt.backend, d)

  # grouping
  groupargs = get(d, :group, nothing) == nothing ? [] : [extractGroupArgs(d[:group], args...)]

  # just in case the backend needs to set up the plot (make it current or something)
  _before_add_series(plt)

  # get the list of dictionaries, one per series
  seriesArgList, xmeta, ymeta = createKWargsList(plt, groupargs..., args...; d...)

  # if we were able to extract guide information from the series inputs, then update the plot
  # @show xmeta, ymeta
  updateDictWithMeta(d, plt.plotargs, xmeta, true)
  updateDictWithMeta(d, plt.plotargs, ymeta, false)

  # now we can plot the series
  for (i,di) in enumerate(seriesArgList)
    plt.n += 1

    if !stringsSupported()
      setTicksFromStringVector(d, di, :x, :xticks)
      setTicksFromStringVector(d, di, :y, :yticks)
    end

    # remove plot args
    for k in keys(_plotDefaults)
      delete!(di, k)
    end

    dumpdict(di, "Series $i")

    _add_series(plt.backend, plt; di...)
  end

  _add_annotations(plt, d)

  warnOnUnsupportedScales(plt.backend, d)

  # add title, axis labels, ticks, etc
  if !haskey(d, :subplot)
    merge!(plt.plotargs, d)
    dumpdict(plt.plotargs, "Updating plot items")
    _update_plot(plt, plt.plotargs)
  end

  _update_plot_pos_size(plt, d)

  current(plt)

  # NOTE: lets ignore the show param and effectively use the semicolon at the end of the REPL statement
  # # do we want to show it?
  if haskey(d, :show) && d[:show]
    gui()
  end

  plt
end

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

# if x or y are a vector of strings, we should create a list of unique strings,
# and map x/y to be the index of the string... then set the x/y tick labels
function setTicksFromStringVector(d::Dict, di::Dict, sym::Symbol, ticksym::Symbol)
  # if the x or y values are strings, set ticks to the unique values, and x/y to the indices of the ticks

  v = di[sym]
  isa(v, AbstractArray) || return

  T = eltype(v)
  if T <: @compat(AbstractString) || (!isempty(T.types) && all(x -> x <: @compat(AbstractString), T.types))

    ticks = unique(di[sym])
    di[sym] = Int[findnext(ticks, v, 1) for v in di[sym]]

    if !haskey(d, ticksym) || d[ticksym] == :auto
      d[ticksym] = (collect(1:length(ticks)), UTF8String[t for t in ticks])
    end
  end
end

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

_before_add_series(plt::Plot) = nothing

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

# should we update the x/y label given the meta info during input slicing?
function updateDictWithMeta(d::Dict, plotargs::Dict, meta::Symbol, isx::Bool)
  lsym = isx ? :xlabel : :ylabel
  if plotargs[lsym] == default(lsym)
    d[lsym] = string(meta)
  end
end
updateDictWithMeta(d::Dict, plotargs::Dict, meta, isx::Bool) = nothing

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

annotations(::@compat(Void)) = []
annotations{X,Y,V}(v::AVec{@compat(Tuple{X,Y,V})}) = v
annotations{X,Y,V}(t::@compat(Tuple{X,Y,V})) = [t]
annotations(anns) = error("Expecting a tuple (or vector of tuples) for annotations: ",
                       "(x, y, annotation)\n    got: $(typeof(anns))")

function _add_annotations(plt::Plot, d::Dict)
  anns = annotations(get(d, :annotation, nothing))
  if !isempty(anns)
    _add_annotations(plt, anns)
  end
end

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

function Base.copy(plt::Plot)
  backend(plt.backend)
  plt2 = plot(; plt.plotargs...)
  for sargs in plt.seriesargs
    sargs = filter((k,v) -> haskey(_seriesDefaults,k), sargs)
    plot!(plt2; sargs...)
  end
  plt2
end

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

# create a new "createKWargsList" which converts all inputs into xs = Any[xitems], ys = Any[yitems].
# Special handling for: no args, xmin/xmax, parametric, dataframes
# Then once inputs have been converted, build the series args, map functions, etc.
# This should cut down on boilerplate code and allow more focused dispatch on type
# note: returns meta information... mainly for use with automatic labeling from DataFrames for now

typealias FuncOrFuncs @compat(Union{Function, AVec{Function}})

# missing
convertToAnyVector(v::@compat(Void); kw...) = Any[nothing], nothing

# fixed number of blank series
convertToAnyVector(n::Integer; kw...) = Any[zeros(0) for i in 1:n], nothing

# numeric vector
convertToAnyVector{T<:Real}(v::AVec{T}; kw...) = Any[v], nothing

# string vector
convertToAnyVector{T<:@compat(AbstractString)}(v::AVec{T}; kw...) = Any[v], nothing

# numeric matrix
convertToAnyVector{T<:Real}(v::AMat{T}; kw...) = Any[v[:,i] for i in 1:size(v,2)], nothing

# function
convertToAnyVector(f::Function; kw...) = Any[f], nothing

# surface
convertToAnyVector(s::Surface; kw...) = Any[s], nothing

# vector of OHLC
convertToAnyVector(v::AVec{OHLC}; kw...) = Any[v], nothing

# dates
convertToAnyVector{D<:Union{Date,DateTime}}(dts::AVec{D}; kw...) = Any[dts], nothing

# list of things (maybe other vectors, functions, or something else)
function convertToAnyVector(v::AVec; kw...)
  if all(x -> typeof(x) <: Real, v)
    # all real numbers wrap the whole vector as one item
    Any[convert(Vector{Float64}, v)], nothing
  else
    # something else... treat each element as an item
    Any[vi for vi in v], nothing
  end
end

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

# in computeXandY, we take in any of the possible items, convert into proper x/y vectors, then return.
# this is also where all the "set x to 1:length(y)" happens, and also where we assert on lengths.
computeX(x::@compat(Void), y) = 1:length(y)
computeX(x, y) = copy(x)
computeY(x, y::Function) = map(y, x)
computeY(x, y) = copy(y)
function computeXandY(x, y)
  if x == nothing && isa(y, Function)
    error("If you want to plot the function `$y`, you need to define the x values somehow!")
  end
  x, y = computeX(x,y), computeY(x,y)
  # @assert length(x) == length(y)
  x, y
end

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

# create n=max(mx,my) series arguments. the shorter list is cycled through
# note: everything should flow through this
function createKWargsList(plt::PlottingObject, x, y; kw...)
  xs, xmeta = convertToAnyVector(x; kw...)
  ys, ymeta = convertToAnyVector(y; kw...)

  mx = length(xs)
  my = length(ys)
  ret = Any[]
  for i in 1:max(mx, my)

    # try to set labels using ymeta
    d = Dict(kw)
    if !haskey(d, :label) && ymeta != nothing
      if isa(ymeta, Symbol)
        d[:label] = string(ymeta)
      elseif isa(ymeta, AVec{Symbol})
        d[:label] = string(ymeta[mod1(i,length(ymeta))])
      end
    end

    # build the series arg dict
    numUncounted = get(d, :numUncounted, 0)
    n = plt.n + i + numUncounted
    dumpdict(d, "before getSeriesArgs")
    d = getSeriesArgs(plt.backend, getplotargs(plt, n), d, i + numUncounted, convertSeriesIndex(plt, n), n)
    dumpdict(d, "after getSeriesArgs")
    d[:x], d[:y] = computeXandY(xs[mod1(i,mx)], ys[mod1(i,my)])

    if haskey(d, :idxfilter)
      d[:x] = d[:x][d[:idxfilter]]
      d[:y] = d[:y][d[:idxfilter]]
    end

    # for linetype `line`, need to sort by x values
    if d[:linetype] == :line
      # order by x
      indices = sortperm(d[:x])
      d[:x] = d[:x][indices]
      d[:y] = d[:y][indices]
      d[:linetype] = :path
    end

    # cleanup those fields that were used only for generating kw args
    for k in (:idxfilter, :numUncounted, :dataframe)
      delete!(d, k)
    end

    # add it to our series list
    push!(ret, d)
  end

  ret, xmeta, ymeta
end

# handle grouping
function createKWargsList(plt::PlottingObject, groupby::GroupBy, args...; kw...)
  ret = Any[]
  for (i,glab) in enumerate(groupby.groupLabels)
    # TODO: don't automatically overwrite labels
    kwlist, xmeta, ymeta = createKWargsList(plt, args...; kw...,
                                            idxfilter = groupby.groupIds[i],
                                            label = string(glab),
                                            numUncounted = length(ret))  # we count the idx from plt.n + numUncounted + i
    append!(ret, kwlist)
  end
  ret, nothing, nothing # TODO: handle passing meta through
end

# pass it off to the x/y version
function createKWargsList(plt::PlottingObject, y; kw...)
  createKWargsList(plt, nothing, y; kw...)
end

# 3d line or scatter
function createKWargsList(plt::PlottingObject, x::AVec, y::AVec, zvec::AVec; kw...)
  d = Dict(kw)
  if !(get(d, :linetype, :none) in _3dTypes)
    d[:linetype] = :path3d
  end
  createKWargsList(plt, x, y; z=zvec, d...)
end

# contours or surfaces... function grid
function createKWargsList(plt::PlottingObject, x::AVec, y::AVec, zf::Function; kw...)
  # only allow sorted x/y for now
  # TODO: auto sort x/y/z properly
  @assert x == sort(x)
  @assert y == sort(y)
  surface = Float64[zf(xi, yi) for xi in x, yi in y]
  createKWargsList(plt, x, y, surface; kw...)  # passes it to the zmat version
end

# contours or surfaces... matrix grid
function createKWargsList{T<:Real}(plt::PlottingObject, x::AVec, y::AVec, zmat::AMat{T}; kw...)
  # only allow sorted x/y for now
  # TODO: auto sort x/y/z properly
  @assert x == sort(x)
  @assert y == sort(y)
  @assert size(zmat) == (length(x), length(y))
  # surf = Surface(convert(Matrix{Float64}, zmat))
  # surf = Array(Any,1,1)
  # surf[1,1] = convert(Matrix{Float64}, zmat)
  d = Dict(kw)
  d[:z] = Surface(convert(Matrix{Float64}, zmat))
  if !(get(d, :linetype, :none) in (:contour, :surface, :wireframe))
    d[:linetype] = :contour
  end
  createKWargsList(plt, x, y; d...) #, z = surf)
end

# contours or surfaces... general x, y grid
function createKWargsList{T<:Real}(plt::PlottingObject, x::AMat{T}, y::AMat{T}, zmat::AMat{T}; kw...)
  @assert size(zmat) == size(x) == size(y)
  surf = Surface(convert(Matrix{Float64}, zmat))
  # surf = Array(Any,1,1)
  # surf[1,1] = convert(Matrix{Float64}, zmat)
  d = Dict(kw)
  d[:z] = Surface(convert(Matrix{Float64}, zmat))
  if !(get(d, :linetype, :none) in (:contour, :surface, :wireframe))
    d[:linetype] = :contour
  end
  createKWargsList(plt, Any[x], Any[y]; d...) #kw..., z = surf, linetype = :contour)
end

function createKWargsList(plt::PlottingObject, surf::Surface; kw...)
  createKWargsList(plt, 1:size(surf.surf,1), 1:size(surf.surf,2), convert(Matrix{Float64}, surf.surf); kw...)
end

function createKWargsList(plt::PlottingObject, x::AVec, y::AVec, surf::Surface; kw...)
  createKWargsList(plt, x, y, convert(Matrix{Float64}, surf.surf); kw...)
end

function createKWargsList(plt::PlottingObject, f::FuncOrFuncs; kw...)
  error("Can't pass a Function or Vector{Function} for y without also passing x")
end

# list of functions
function createKWargsList(plt::PlottingObject, f::FuncOrFuncs, x; kw...)
  @assert !(typeof(x) <: FuncOrFuncs)  # otherwise we'd hit infinite recursion here
  createKWargsList(plt, x, f; kw...)
end

# special handling... xmin/xmax with function(s)
function createKWargsList(plt::PlottingObject, f::FuncOrFuncs, xmin::Real, xmax::Real; kw...)
  width = get(plt.plotargs, :size, (100,))[1]
  x = collect(linspace(xmin, xmax, width))  # we don't need more than the width
  createKWargsList(plt, x, f; kw...)
end

mapFuncOrFuncs(f::Function, u::AVec) = map(f, u)
mapFuncOrFuncs(fs::AVec{Function}, u::AVec) = [map(f, u) for f in fs]

# special handling... xmin/xmax with parametric function(s)
createKWargsList{T<:Real}(plt::PlottingObject, fx::FuncOrFuncs, fy::FuncOrFuncs, u::AVec{T}; kw...) = createKWargsList(plt, mapFuncOrFuncs(fx, u), mapFuncOrFuncs(fy, u); kw...)
createKWargsList{T<:Real}(plt::PlottingObject, u::AVec{T}, fx::FuncOrFuncs, fy::FuncOrFuncs; kw...) = createKWargsList(plt, mapFuncOrFuncs(fx, u), mapFuncOrFuncs(fy, u); kw...)
createKWargsList(plt::PlottingObject, fx::FuncOrFuncs, fy::FuncOrFuncs, umin::Real, umax::Real, numPoints::Int = 1000; kw...) = createKWargsList(plt, fx, fy, linspace(umin, umax, numPoints); kw...)

# special handling... no args... 1 series
function createKWargsList(plt::PlottingObject; kw...)
  d = Dict(kw)
  if !haskey(d, :y)
    # assume we just want to create an empty plot object which can be added to later
    return [], nothing, nothing
    # error("Called plot/subplot without args... must set y in the keyword args.  Example: plot(; y=rand(10))")
  end

  if haskey(d, :x)
    return createKWargsList(plt, d[:x], d[:y]; kw...)
  else
    return createKWargsList(plt, d[:y]; kw...)
  end
end

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

"For DataFrame support.  Imports DataFrames and defines the necessary methods which support them."
function dataframes()
  @eval import DataFrames

  @eval function createKWargsList(plt::PlottingObject, df::DataFrames.AbstractDataFrame, args...; kw...)
    createKWargsList(plt, args...; kw..., dataframe = df)
  end

  # expecting the column name of a dataframe that was passed in... anything else should error
  @eval function extractGroupArgs(s::Symbol, df::DataFrames.AbstractDataFrame, args...)
    if haskey(df, s)
      return extractGroupArgs(df[s])
    else
      error("Got a symbol, and expected that to be a key in d[:dataframe]. s=$s d=$d")
    end
  end

  @eval function getDataFrameFromKW(; kw...)
    for (k,v) in kw
      if k == :dataframe
        return v
      end
    end
    error("Missing dataframe argument in arguments!")
  end

  # the conversion functions for when we pass symbols or vectors of symbols to reference dataframes
  @eval convertToAnyVector(s::Symbol; kw...) = Any[getDataFrameFromKW(;kw...)[s]], s
  @eval convertToAnyVector(v::AVec{Symbol}; kw...) = (df = getDataFrameFromKW(;kw...); Any[df[s] for s in v]), v
end

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


text:
#######################################################################################################################################################################################################
#
# Changes to this function
# General
#     2021-Aug-04: add this function to acoid NaN calculations
#
#######################################################################################################################################################################################################
"""
    nanmean(x::Array)

Return the mean of array by ommiting the NaN, given
- `x` Array of numbers, can be NaN

---
Example
```julia
xs = [1, 2, 4, NaN];
nmean = nanmean(xs);
```
"""
function nanmean(x::Array)
    _x = filter(!isnan, x);

    if length(_x) == 0 return NaN end;

    return mean( _x )
end


text:
# -*- coding: utf-8 -*-
# ---
# jupyter:
#   jupytext:
#     formats: ipynb,jl:hydrogen
#     text_representation:
#       extension: .jl
#       format_name: hydrogen
#       format_version: '1.3'
#       jupytext_version: 1.11.2
#   kernelspec:
#     display_name: Julia 1.8.0-DEV
#     language: julia
#     name: julia-1.8
# ---

# %% [markdown]
# https://twitter.com/yujitach/status/1424030835771023363

# %%
using LinearAlgebra
using LoopVectorization
using LinearMaps
using ArnoldiMethod

function prepareDiag(L)
    diag = zeros(2^L)
    for state = 1:2^L
        for i = 1:L
            j = i==L ? 1 : i+1
            @inbounds diag[state] -= (((state >> (i-1))&1) == ((state >> (j-1))&1)) ? 1 : -1
        end
    end
    diag
end

function Hfunc!(C, B, diag, L)
    N = length(diag)
    @tturbo for state = 1:N
        C[state] = diag[state] * B[state]
    end
    for i = 1:L
        @tturbo for state = 1:N
            newstate = (state&(~(2^L))) ⊻ (1<<(i-1))
            c = newstate == 0
            newstate = !c*newstate + c*N # remove if statement
            C[newstate] -= B[state]
        end
    end
end
prepareHfunc!(L) = (C, B) -> Hfunc!(C, B, prepareDiag(L), L)

function f_LinearMap(L)
    H_LinearMap = LinearMap(prepareHfunc!(L), 2^L, ismutating=true, issymmetric=true, isposdef=false)
    d_LM, h_LM = partialschur(H_LinearMap, nev=1, which=SR())
    d_LM.eigenvalues[1]
end

f_exact(L) = -2sum(abs(sin((n-1/2) * pi/L)) for n in 1:L)

f_LinearMap(20), f_exact(20)

# %%
@time f_LinearMap(20)
@time f_LinearMap(20)
@time f_LinearMap(20)

# %% [markdown]
# https://juliaphysics.github.io/PhysicsTutorials.jl/tutorials/general/quantum_ising/quantum_ising.html

# %%
using LinearAlgebra
using SparseArrays
using ArnoldiMethod
⊗(x,y) = kron(x,y)

function TransverseFieldIsing_sparse(;N,h)
    id = [1 0; 0 1] |> sparse
    σˣ = [0 1; 1 0] |> sparse
    σᶻ = [1 0; 0 -1] |> sparse

    first_term_ops = fill(id, N)
    first_term_ops[1] = σᶻ
    first_term_ops[2] = σᶻ

    second_term_ops = fill(id, N)
    second_term_ops[1] = σˣ

    H = spzeros(Int, 2^N, 2^N)
    for i in 1:N
        H -= foldl(⊗, first_term_ops)
        first_term_ops = circshift(first_term_ops,1)
    end

    for i in 1:N
        H -= h*foldl(⊗, second_term_ops)
        second_term_ops = circshift(second_term_ops,1)
    end
    H
end

function f_sparse(L)
    H = TransverseFieldIsing_sparse(N=L, h=1)
    d, hist = partialschur(H; nev=1, which=SR())
    d.eigenvalues[1]
end

f_exact(L) = -2sum(abs(sin((n-1/2) * pi/L)) for n in 1:L)

f_sparse(20), f_exact(20)

# %%
L = 20
@time H = TransverseFieldIsing_sparse(N=L, h=1)
@time d, hist = partialschur(H; nev=1, which=SR())

# %%
L = 20
@time H = TransverseFieldIsing_sparse(N=L, h=1)
@time d, hist = partialschur(H; nev=1, which=SR())

# %%
L = 20
@time H = TransverseFieldIsing_sparse(N=L, h=1)
@time d, hist = partialschur(H; nev=1, which=SR())

# %%
TransverseFieldIsing_sparse(N=20, h=1)

# %%
function TransverseFieldIsing_sparse_revised(;N, h)
    id = [1 0; 0 1] |> sparse
    σˣ = [0 1; 1 0] |> sparse
    σᶻ = [1 0; 0 -1] |> sparse

    first_term_ops = fill(id, N)
    first_term_ops[1] = σᶻ
    first_term_ops[2] = σᶻ

    second_term_ops = fill(id, N)
    second_term_ops[1] = σˣ

    tmp = [spzeros(Int, 2^k, 2^k) for k in 1:N]
    H = spzeros(Int, 2^N, 2^N)

    for i in 1:N
        tmp[1] .= first_term_ops[mod1(i, N)]
        for k in 2:N
            kron!(tmp[k], tmp[k-1], first_term_ops[mod1(i+k-1, N)])
        end
        H .-= tmp[N]
    end

    for i in 1:N
        tmp[1] .= second_term_ops[mod1(i, N)]
        for k in 2:N
            kron!(tmp[k], tmp[k-1], second_term_ops[mod1(i+k-1, N)])
        end
        H .-= h .* tmp[N]
    end
    H
end

function f_sparse_revised(L)
    H = TransverseFieldIsing_sparse_revised(N=L, h=1)
    d, hist = partialschur(H; nev=1, which=SR())
    d.eigenvalues[1]
end

TransverseFieldIsing_sparse_revised(N=10, h=1) == TransverseFieldIsing_sparse(N=10, h=1)

# %%
L = 20
@time H = TransverseFieldIsing_sparse_revised(N=L, h=1)
@time d, hist = partialschur(H; nev=1, which=SR())

# %%
L = 20
@time H = TransverseFieldIsing_sparse_revised(N=L, h=1)
@time d, hist = partialschur(H; nev=1, which=SR())

# %%
L = 20
@time H = TransverseFieldIsing_sparse_revised(N=L, h=1)
@time d, hist = partialschur(H; nev=1, which=SR())

# %%
TransverseFieldIsing_sparse_revised(N=20, h=1)

# %%
bit(a, k) = (a >> k) & 1

function sigmaz2(N)
    d = zeros(Int, 2^N)
    for k in 1:N, a in 0:2^N-1
        d[a+1] -= ifelse(bit(a, N-k) == bit(a, N - mod1(k+1, N)), 1, -1)
    end
    d
end

function sigmax(N, h, k)
    v = zeros(Int, 2^N - 2^(N-k))
    for i in 1:2:2^k-1
        v[(2^(N-k)*(i-1) + 1):2^(N-k)*i] .= -h
    end
    v
end

function TransverseFieldIsing_sparse_revised2(; N, h)
    d = sigmaz2(N)
    v = sigmax.(N, h, 1:N)
    H = spdiagm(
        (-2^(N-k) => v[k] for k in 1:N)...,
        0 => d,
        ( 2^(N-k) => v[k] for k in 1:N)...)
end

function f_sparse_revised2(L)
    H = TransverseFieldIsing_sparse_revised2(N=L, h=1)
    d, hist = partialschur(H; nev=1, which=SR())
    d.eigenvalues[1]
end

TransverseFieldIsing_sparse_revised2(N=10, h=1) == TransverseFieldIsing_sparse(N=10, h=1)

# %%
L = 20
@time H = TransverseFieldIsing_sparse_revised2(N=L, h=1)
@time d, hist = partialschur(H; nev=1, which=SR())

# %%
L = 20
@time H = TransverseFieldIsing_sparse_revised2(N=L, h=1)
@time d, hist = partialschur(H; nev=1, which=SR())

# %%
L = 20
@time H = TransverseFieldIsing_sparse_revised2(N=L, h=1)
@time d, hist = partialschur(H; nev=1, which=SR())

# %%
function TransverseFieldIsing_sparse_revised2_time(;N, h)
    @time d = sigmaz2(N)
    @time v = sigmax.(N, h, 1:N)
    @time H = spdiagm(
        (-2^(N-k) => v[k] for k in 1:N)...,
        0 => d,
        ( 2^(N-k) => v[k] for k in 1:N)...)
end

TransverseFieldIsing_sparse_revised2_time(N=20, h=1); println()
TransverseFieldIsing_sparse_revised2_time(N=20, h=1); println()
TransverseFieldIsing_sparse_revised2_time(N=20, h=1)

# %%
bit(a, k) = (a >> k) & 1

function sigmaz2_sparse(N)
    d = zeros(Int, 2^N)
    for k in 1:N, a in 0:2^N-1
        d[a+1] -= ifelse(bit(a, N-k) == bit(a, N - mod1(k+1, N)), 1, -1)
    end
    sparse(d) # should be sparse
end

function sigmax_sparse(N, h, k)
    v = zeros(Int, 2^N - 2^(N-k))
    for i in 1:2:2^k-1
        v[(2^(N-k)*(i-1) + 1):2^(N-k)*i] .= -h
    end
    sparse(v) # should be sparse
end

function TransverseFieldIsing_sparse_revised3(; N, h)
    d = sigmaz2_sparse(N)
    v = sigmax_sparse.(N, h, 1:N)
    H = spdiagm(
        (-2^(N-k) => v[k] for k in 1:N)...,
        0 => d,
        ( 2^(N-k) => v[k] for k in 1:N)...)
end

function f_sparse_revised3(L)
    H = TransverseFieldIsing_sparse_revised3(N=L, h=1)
    d, hist = partialschur(H; nev=1, which=SR())
    d.eigenvalues[1]
end

TransverseFieldIsing_sparse_revised3(N=10, h=1) == TransverseFieldIsing_sparse(N=10, h=1)

# %% tags=[]
L = 20
@time H = TransverseFieldIsing_sparse_revised3(N=L, h=1)
@time d, hist = partialschur(H; nev=1, which=SR())

# %%
L = 20
@time H = TransverseFieldIsing_sparse_revised3(N=L, h=1)
@time d, hist = partialschur(H; nev=1, which=SR())

# %%
L = 20
@time H = TransverseFieldIsing_sparse_revised3(N=L, h=1)
@time d, hist = partialschur(H; nev=1, which=SR())

# %%
function TransverseFieldIsing_sparse_revised3_time(; N, h)
    @time d = sigmaz2_sparse(N)
    @time v = sigmax_sparse.(N, h, 1:N)
    @time H = spdiagm(
        (-2^(N-k) => v[k] for k in 1:N)...,
        0 => d,
        ( 2^(N-k) => v[k] for k in 1:N)...)
end

TransverseFieldIsing_sparse_revised3_time(N=20, h=1); println()
TransverseFieldIsing_sparse_revised3_time(N=20, h=1); println()
TransverseFieldIsing_sparse_revised3_time(N=20, h=1)

# %%
bit(a, k) = (a >> k) & 1

function sigmaz2_sparse(N)
    d = zeros(Int8, 2^N)
    for k in 1:N, a in 0:2^N-1
        d[a+1] -= ifelse(bit(a, N-k) == bit(a, N - mod1(k+1, N)), 1, -1)
    end
    sparse(d)
end

function IJV(N, h=1)
    i, v = findnz(sigmaz2_sparse(N))
    m = length(i)
    M = 2^N

    I = Vector{Int32}(undef, m + N*M)
    J = similar(I)
    V = similar(I, Int8)

    I[1:m] .= i
    J[1:m] .= i
    V[1:m] .= v
    V[m+1:end] .= -h

    for a in 0:M-1
        s = m + N*a
        I[s+1:s+N] .= a+1
        for i in 1:N
            b = a ⊻ (1 << (N-i))
            J[s+i] = b+1
        end
    end

    I, J, V
end

TransverseFieldIsing_sparse_revised4(; N, h) = sparse(IJV(N, h)...)

function f_sparse_revised4(L)
    H = TransverseFieldIsing_sparse_revised4(N=L, h=1)
    d, hist = partialschur(H; nev=1, which=SR())
    d.eigenvalues[1]
end

TransverseFieldIsing_sparse_revised4(N=10, h=1) == TransverseFieldIsing_sparse(N=10, h=1)

# %%
L = 20
@time H = TransverseFieldIsing_sparse_revised4(N=L, h=1)
@time d, hist = partialschur(H; nev=1, which=SR())

# %%
L = 20
@time H = TransverseFieldIsing_sparse_revised4(N=L, h=1)
@time d, hist = partialschur(H; nev=1, which=SR())

# %%
L = 20
@time H = TransverseFieldIsing_sparse_revised4(N=L, h=1)
@time d, hist = partialschur(H; nev=1, which=SR())

# %%
function TransverseFieldIsing_sparse_revised4_time(; N, h)
    @time I, J, V = IJV(N, h)
    @time sparse(I, J, V)
end

TransverseFieldIsing_sparse_revised4_time(N=20, h=1); println()
TransverseFieldIsing_sparse_revised4_time(N=20, h=1); println()
TransverseFieldIsing_sparse_revised4_time(N=20, h=1)

# %%

max_stars_repo_path: 0015/Arpack/Sparse version revised.jl

max_stars_repo_name: genkuroki/public

max_issues_repo_path: 0015/Arpack/Sparse version revised.jl

max_issues_repo_name: genkuroki/public

max_forks_repo_path: 0015/Arpack/Sparse version revised.jl

max_forks_repo_name: genkuroki/public


text:
using Treesitter

p = Parser(tree_sitter_c())
tree = parse_string(p, "typedef int* jeff;")

sexp = repr("text/plain",tree.root_node)

l1 = node_children(tree.root_node)
l2 = node_children.(l1)


text:
dictkeys(d::Dict) = (collect(keys(d))...,)
dictvalues(d::Dict) = (collect(values(d))...,)

#=
to_nt(d::Dict{Symbol,T}) where {T} =
    NamedTuple{dictkeys(d)}(dictvalues(d))
=#
typed(v::Number) = Float64(v)
function typed(v::String)
  v_float = tryparse(Float64,v)
  ret = isnothing(v_float) ? Symbol(v) : v_float
  ret
end

typed(v::Symbol) = v

function _subset(
  vector::AbstractVector,
  subset::AbstractVector{P} = Pair{Symbol, Symbol}[]
) where P <: Pair{Symbol, Symbol}
  if length(subset) === 0
    # use the whole vector
    return vector
  else
    # this function performs matching
    # returs true if all values in subset equal to values in row
    subset_fun = (row) -> begin
      res = [row[key] == value for (key, value) in subset]
      return all(res)
    end

    return filter(subset_fun, vector)
  end
end

has_saveat(c::Scenario) = !isempty(c.prob.kwargs[:callback].discrete_callbacks[1].affect!.saveat_cache)
has_saveat(mcr::MCResult) = mcr.saveat

# auxilary method to update values in LVector
# TODO: not optimal method
# TODO: write for any number of pairs update(la::LArray, argv...)
function update(la::LArray, pairs::AbstractVector{Pair{Symbol, C}}) where C<:Real
  clone = copy(la)
  keys_la = keys(la)

  for (key, value) in pairs
    key ∉ keys_la && throw("key :$key is not found in LVector.")
    clone[key] = value
  end

  return clone
end

update(p1::AbstractVector{P}, p2::AbstractVector{P}) where P<:Pair = update(Dict(p1),Dict(p2))

function update(d1::Dict, d2::Dict)
  d = copy(d1)
  keys_d = keys(d)

  for (k,v) in d2
    k ∉ keys_d && throw("key :$k is not found in Dict.")
    d[k] = v
  end

  return d
end

float64(n::Number) = Float64(n)
float64(m::Missing) = m

bool(i::Int) = Bool(i)
function bool(s::AbstractString)
  s == "TRUE" && return true
  s == "FALSE" && return false
  return parse(Bool, s)
end
bool(b::Bool) = b