text:
#
# Code generation for intrinsic functions.
#

make_intrinsic_arg!(cg, x) = codegen!(cg, x)
make_intrinsic_arg!(cg, x::Type{T}) where {T} = llvmtype(x)
make_intrinsic_arg!(cg, x::GlobalRef) = make_intrinsic_arg!(cg, eval(x))

function uint_cnvt!(cg, to, x)
    t = LLVM.llvmtype(x)
    t == to && return x
    width(to) < width(t) && return LLVM.trunc!(cg.builder, x, to)
    return LLVM.zext!(cg.builder, x, to)
end

function emit_intrinsic!(cg::CodeCtx, name, jlargs)
    args = Any[]
    for v in jlargs
        push!(args, make_intrinsic_arg!(cg, v))
    end
    name == :neg_int  && return LLVM.neg!(cg.builder, args[1])
    name == :add_int  && return LLVM.add!(cg.builder, args[1], args[2])
    name == :sub_int  && return LLVM.sub!(cg.builder, args[1], args[2])
    name == :mul_int  && return LLVM.mul!(cg.builder, args[1], args[2])
    name == :sdiv_int && return LLVM.sdiv!(cg.builder, args[1], args[2])
    name == :udiv_int && return LLVM.udiv!(cg.builder, args[1], args[2])
    name == :srem_int && return LLVM.srem!(cg.builder, args[1], args[2])
    name == :urem_int && return LLVM.urem!(cg.builder, args[1], args[2])
    name == :neg_float  && return LLVM.fsub!(cg.builder, codegen!(cg, -0.0), args[1])
    name == :add_float  && return LLVM.fadd!(cg.builder, args[1], args[2])
    name == :sub_float  && return LLVM.fsub!(cg.builder, args[1], args[2])
    name == :mul_float  && return LLVM.fmul!(cg.builder, args[1], args[2])
    name == :div_float  && return LLVM.fdiv!(cg.builder, args[1], args[2])
    name == :rem_float  && return LLVM.frem!(cg.builder, args[1], args[2])
    ## WRONG. Not sure how to do these. Need a new builder with fast math set.
    name == :neg_float_fast  && return LLVM.fneg!(cg.builder, args[1])
    name == :add_float_fast  && return LLVM.fadd!(cg.builder, args[1], args[2])
    name == :sub_float_fast  && return LLVM.fsub!(cg.builder, args[1], args[2])
    name == :mul_float_fast  && return LLVM.fmul!(cg.builder, args[1], args[2])
    name == :div_float_fast  && return LLVM.fdiv!(cg.builder, args[1], args[2])
    name == :rem_float_fast  && return LLVM.frem!(cg.builder, args[1], args[2])
    ## More tough ones
    # name == :fma_float  && return LLVM.frem!(cg.builder, args[1], args[2], "sremtmp")

    function ltyp(x)
        res = string(LLVM.llvmtype(x))
        if res == "double"
            res = "f64"
        elseif res == "float"
            res = "f32"
        end
        return res
    end

    name == :muladd_float  && return LLVM.call!(cg.builder,
        cg.extern[Symbol("llvm.fmuladd.$(ltyp(args[1]))")],
        LLVM.Value[args[1], args[2], args[3]])

    function checked_int(sym)
        ifun = cg.extern[Symbol("llvm.$sym.with.overflow.$(ltyp(args[1]))")]
        res = LLVM.call!(cg.builder, ifun, LLVM.Value[args[1], args[2]])
        val = LLVM.extract_value!(cg.builder, res, 0)
        obit = LLVM.extract_value!(cg.builder, res, 1)
        obyte = LLVM.zext!(cg.builder, obit, int8_t)
        tupptr = alloca!(cg.builder, LLVM.StructType([LLVM.llvmtype(args[1]), int8_t], CodeGen.ctx))
        store!(cg, val,   LLVM.struct_gep!(cg.builder, tupptr, 0))
        store!(cg, obyte, LLVM.struct_gep!(cg.builder, tupptr, 1))
        return LLVM.load!(cg.builder, tupptr)
    end

    name == :checked_sadd_int  && return checked_int(:sadd)
    name == :checked_uadd_int  && return checked_int(:uadd)
    name == :checked_ssub_int  && return checked_int(:ssub)
    name == :checked_usub_int  && return checked_int(:usub)
    name == :checked_smul_int  && return checked_int(:smul)
    name == :checked_umul_int  && return checked_int(:umul)
    # name == :checked_sdiv_int  && return LLVM.sdiv!(cg.builder, args[1], args[2])
    # name == :checked_udiv_int  && return LLVM.udiv!(cg.builder, args[1], args[2])
    # name == :checked_srem_int  && return LLVM.srem!(cg.builder, args[1], args[2])
    # name == :checked_urem_int  && return LLVM.urem!(cg.builder, args[1], args[2])
    ##

    emit_bool!(x) = LLVM.zext!(cg.builder, x, int8_t)
    emit_icmp!(fun) = emit_bool!(LLVM.icmp!(cg.builder, fun, args[1], args[2]))
    emit_fcmp!(fun) = emit_bool!(LLVM.fcmp!(cg.builder, fun, args[1], args[2]))
    name == :eq_int  && return emit_icmp!(LLVM.API.LLVMIntEQ)
    name == :ne_int  && return emit_icmp!(LLVM.API.LLVMIntNE)
    name == :slt_int && return emit_icmp!(LLVM.API.LLVMIntSLT)
    name == :ult_int && return emit_icmp!(LLVM.API.LLVMIntULT)
    name == :sle_int && return emit_icmp!(LLVM.API.LLVMIntSLE)
    name == :ule_int && return emit_icmp!(LLVM.API.LLVMIntULE)
    name == :eq_float && return emit_fcmp!(LLVM.API.LLVMRealOEQ)
    name == :ne_float && return emit_fcmp!(LLVM.API.LLVMRealONE)
    name == :lt_float && return emit_fcmp!(LLVM.API.LLVMRealOLT)
    name == :le_float && return emit_fcmp!(LLVM.API.LLVMRealOLE)
    ## need fast versions of above
    # fpiseq
    # fpislt
    name == :and_int  && return emit_bool!(LLVM.and!(cg.builder, args[1], args[2]))
    name == :or_int   && return emit_bool!(LLVM.or!(cg.builder, args[1], args[2]))
    name == :xor_int  && return emit_bool!(LLVM.xor!(cg.builder, args[1], args[2]))

    a1(name) = LLVM.call!(cg.builder, cg.extern[Symbol("llvm.$name.$(ltyp(args[1]))")], LLVM.Value[args[1]])

    name == :bswap_int  && return a1("bswap")
    name == :ctpop_int  && return a1("ctpop")
    name == :abs_float  && return a1("fabs")
    name == :ceil_llvm  && return a1("ceil")
    name == :floor_llvm && return a1("floor")
    name == :trunc_llvm && return a1("trunc")
    name == :rint_llvm  && return a1("rint")
    name == :sqrt_llvm  && return a1("sqrt")

    name == :ctlz_int  && return LLVM.call!(cg.builder,
        cg.extern[Symbol("llvm.ctlz.$(ltyp(args[1]))")],
        LLVM.Value[args[1], LLVM.ConstantInt(int1_t, 0)])

    # cttz_int
    # copysign_float
    if name == :flipsign_int
        # ignores the constant case
        tmp = LLVM.ashr!(cg.builder, args[2], codegen!(cg, LLVM.width(LLVM.llvmtype(args[1])) - 1))
        return LLVM.xor!(cg.builder, LLVM.add!(cg.builder, args[1], tmp), tmp)
    end

    name == :sitofp  && return LLVM.sitofp!(cg.builder, args[2], args[1])
    name == :uitofp  && return LLVM.uitofp!(cg.builder, args[2], args[1])
    name == :fptosi  && return LLVM.fptosi!(cg.builder, args[2], args[1])
    name == :fptoui  && return LLVM.fptoui!(cg.builder, args[2], args[1])
    name == :trunc_int  && return LLVM.trunc!(cg.builder, args[2], args[1])
    name == :fptrunc    && return LLVM.fptrunc!(cg.builder, args[2], args[1])
    name == :fpext      && return LLVM.fpext!(cg.builder, args[2], args[1])
    name == :sext_int   && return LLVM.sext!(cg.builder, args[2], args[1])
    name == :zext_int   && return LLVM.zext!(cg.builder, args[2], args[1])
    name == :fpzext     && return LLVM.fpext!(cg.builder, args[2], args[1])
    if name == :arraylen
        p = LLVM.bitcast!(cg.builder, args[1], LLVM.PointerType(llvmtype(Tuple{Ptr{Nothing}, Csize_t})))
        len_p = LLVM.struct_gep!(cg.builder, p, 1)
        return LLVM.load!(cg.builder, len_p)
    end

    function emit_sh!(fun)
        t1 = LLVM.llvmtype(args[1])
        t2 = LLVM.llvmtype(args[2])
        return LLVM.select!(cg.builder,
                            LLVM.icmp!(cg.builder, LLVM.API.LLVMIntUGE, args[2], LLVM.ConstantInt(t2, width(t1))),
                            LLVM.ConstantInt(t1, 0),
                            fun(cg.builder, args[1], uint_cnvt!(cg, t1, args[2])))
    end
    name == :shl_int  && return emit_sh!(LLVM.shl!)
    name == :ashr_int && return emit_sh!(LLVM.ashr!)
    name == :lshr_int && return emit_sh!(LLVM.lshr!)

    ######
    ## NOT UNIVERSAL, but maybe useful for getting a few things to work
    name == :select_value && return LLVM.select!(cg.builder, emit_condition!(cg, args[1]), args[2], args[3])
    name == :not_int      && return emit_bool!(LLVM.not!(cg.builder, emit_condition!(cg, args[1])))
    name == :bitcast      && return LLVM.bitcast!(cg.builder, args[2], args[1])  # not completely general

    ######
    name == :add_ptr && return LLVM.ptrtoint!(cg.builder,
                                              LLVM.gep!(cg.builder,
                                                        LLVM.inttoptr!(cg.builder, args[1], int8_t_ptr),
                                                        [args[2]]),
                                              LLVM.llvmtype(args[1]))
    name == :sub_ptr && return LLVM.ptrtoint!(cg.builder,
                                              LLVM.gep!(cg.builder,
                                                        LLVM.inttoptr!(cg.builder, args[1], int8_t_ptr),
                                                        [LLVM.neg!(cg.builder, args[2])]),
                                              LLVM.llvmtype(args[1]))

    ######
    # name == :pointerref && return LLVM.call!(cg.builder, cg.extern[:jl_pointerref], LLVM.Value[emit_box!(cg, args[1]), emit_box!(cg, args[2]), emit_box!(cg, args[3])])
    # name == :pointerset && return LLVM.call!(cg.builder, cg.extern[:jl_pointerset], LLVM.Value[emit_box!(cg, args[1]), emit_box!(cg, args[2]), emit_box!(cg, args[3]), emit_box!(cg, args[4])])
    name == :pointerref && return LLVM.call!(cg.builder, cg.extern[:jl_pointerref], LLVM.Value[emit_box!(cg, args[1]), emit_box!(cg, args[2]), emit_box!(cg, args[3])])
    name == :pointerset && return LLVM.call!(cg.builder, cg.extern[:jl_pointerset], LLVM.Value[emit_box!(cg, args[1]), emit_box!(cg, args[2]), emit_box!(cg, args[3]), emit_box!(cg, args[4])])

    error("Unsupported intrinsic: $name")
end


text:
using Plots, UnPack, FFTW
using MOONS.Grids
using MOONS.Fields
using MOONS.FFTs

FT = Float64
n = 128
cx = Coordinates(FT(0),FT(1),n)
cy = Coordinates(FT(0),FT(1),n)
cz = Coordinates(FT(0),FT(1),2)
g = Grid((cx,cy,cz))

f_exact = CellCenter(g)
f_recreate = CellCenter(g)
u_exact = CellCenter(g)
f = CellCenter(g)
u = CellCenter(g)

@inbounds for local_grid in GridField(g,f_exact)
    @unpack ijk,x,y = local_grid
    u_exact[ijk] = cos(3*π*x)*cos(3*π*y)
    f_exact[ijk] = -2*(3*π)^2*cos(3*π*x)*cos(3*π*y)
end
∇²!(f, u_exact, g)
extrap!(f, g)
solve!(u, f_exact, g, 3)

∇²!(f_recreate, u, g)

@unpack xc, yc = all_coords(g)

contourf(xc,yc, midslice(f, 3); c = :viridis)
savefig("f_discrete.svg")

contourf(xc,yc, midslice(f_exact, 3); c = :viridis)
savefig("f_exact.svg")

u_ave = average(u, 3);
@show max(u_ave...)
max(abs.(diff(u; dims=3))...) ≈ 0 || @warn "solution is not 2D"

contourf(xc,yc, u_ave; c = :viridis)
savefig("u_discrete.svg")

contourf(xc,yc, midslice(u_exact, 3); c = :viridis, clims=(-1,1))
savefig("u_exact.svg")

u_err = CellCenter(g)
u_err .= abs.(u_exact - u);
assign!(u_err, 0, Ghost())
@show max(u_err...)

contourf(xc,yc, midslice(u_err, 3); c = :viridis)
savefig("u_err.svg")

ext: jl


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

using Markdown
using InteractiveUtils

# ╔═╡ 9f18ccd0-8e25-11eb-20e6-ab941f663550
using Plots

# ╔═╡ 10a9374c-8e1a-11eb-10b4-fb8272dd46d1
md"""
## SIR Model Simulation

#### Intro

The **SIR Model** has been developed to simulate an epidemic over time. The model consists of a system of 3 differential equations that express the rates of change of 3 variables over time. The 3 variables are:

1. **S** - the susceptibles of getting the infection
2. **I** - the infected
3. **R** - the recovered from the infection
3. D - the number of dead

#### The Model
Here follow the 3 equations that govern the model dynamics:

${ \begin{align*}
\frac{\mathrm{d}S}{\mathrm{d}t} & = -\beta \cdot S \cdot I + \delta \cdot R \\
\frac{\mathrm{d}I}{\mathrm{d}t} & = \beta \cdot S \cdot I + \gamma \cdot I \\
\frac{\mathrm{d}R}{\mathrm{d}t} & = \gamma \cdot I - \delta \cdot R \\
\end{align*} }$

#### Model simulation
The following is a simulation of the model described above
"""

# ╔═╡ 43970fa0-8e1b-11eb-3052-1701966a4478
begin
	#Model parameters
	beta1 = 0.14 	#rate of infection
	gamma1 = 0.07 	#rate of recovery (try also 0.07)
	delta1 = 0.00 	#rate of immunity loss
	rho1 = 0.001 	#rate of death of infected
	N1 = 6*10^7 	#Total population N = S + I + R + D
	I0 = 100.0 	#initial number of infected
	beta2 = 0.0 	#rate of infection
	gamma2 = 0.035 	#rate of recovery (try also o.07)
	delta2 = 0.0 	#rate of immunity loss
	rho2 = 0.2 	#rate of death of infected
	beta3 = 0.0 	#rate of infection from group 1 to group 2
	beta4 = 0.0 	#rate of infection from group 2 to group 1
	N2 = 1*10^6 #Total population N = S + I + R + D
	J0 = 0.0 	#initial number of infected
	gamma = 0.0
	T = 38*30 	#period of 300 daysplot
	dt = 1/400 	#time interval of 6 hours (1/4 of a day)

	md"""
	Value of parameter R10 is $(beta1/gamma1) | Value of parameter R20 is $(beta2/gamma2)
	"""
end

# ╔═╡ 883a4ad8-8e26-11eb-06c4-e30f3c0f70c1
# Plots that display the epidemic outbreak
tt = 0:dt:T-dt;

# ╔═╡ 5f948160-8e27-11eb-1003-fbcf284800e6
#Map
#plot(tt,[I1,I2,D1,D2],color=[:blue :red :green :yellow], lw=2, grid=true, xlabel="Days", ylabel="Number of individuals", label=["I1" "I2" "D1" "D2"])

# ╔═╡ c30b3fb4-8e1c-11eb-05c7-d7f95aa97fcf
md"Function that calculates the model evolution at time period T"

# ╔═╡ 01a236b0-8e1d-11eb-2630-ebff3f351bb3
function sir_model(beta1,gamma1,delta1,rho1,N1,beta2,gamma2,delta2,rho2,N2,beta3,beta4,I0,J0,T,dt)
	S1 = zeros(Int(T/dt))
	S1[1] = N1
	I1 = zeros(Int(T/dt))
	I1[1] = I0
	R1 = zeros(Int(T/dt))
	D1 = zeros(Int(T/dt))
	S2 = zeros(Int(T/dt))
	S2[1] = N2
	I2 = zeros(Int(T/dt))
	I2[1] = J0
	R2 = zeros(Int(T/dt))
	D2 = zeros(Int(T/dt))
	I = zeros(Int(T/dt))
	J1 = zeros(Int(T/dt))
	J2 = zeros(Int(T/dt))
	gamma = gamma1

	for tt ∈ 1:Int(T/dt)-1
		if I2[tt] + I1[tt] > 10
			if tt > 10/dt
				if J2[tt-Int(5/dt)] > 200
					gamma = 2.8*gamma1
				end
				if J2[tt-Int(5/dt)] < 17
					gamma = gamma1
				end
			end

			dS1 = (-beta1*I1[tt]*S1[tt]/N1+delta1*R1[tt]-beta4*S1[tt]*I2[tt]/N1) * dt
			dI1 = (beta1*I1[tt]*S1[tt]/N1-gamma*I1[tt]+beta4*S1[tt]*I2[tt]/N1) * dt
			dR1 = (gamma*(1-rho1)*I1[tt]-delta1*R1[tt]) * dt
			dD1 = (gamma*rho1*I1[tt]) * dt
			dJ1 = beta1*I1[tt]*S1[tt]/N1

			S1[tt+1] = S1[tt] + dS1
			I1[tt+1] = I1[tt] + dI1
			R1[tt+1] = R1[tt] + dR1
			D1[tt+1] = D1[tt] + dD1
			J1[tt+1] = J1[tt] + dJ1

			dS2 = (-beta2*I2[tt]*S2[tt]/N2+delta2*R2[tt]-beta3*I1[tt]*S2[tt]/N1) * dt
			dI2 = (beta2*I2[tt]*S2[tt]/N2-gamma2*I2[tt]+beta3*I1[tt]*S2[tt]/N1) * dt
			dR2 = (gamma2*(1-rho2)*I2[tt] - delta2*R2[tt]) * dt
			dD2 = (gamma2*rho2*I2[tt]) * dt

			S2[tt+1] = S2[tt] + dS2
			I2[tt+1] = I2[tt] + dI2
			R2[tt+1] = R2[tt] + dR2
			D2[tt+1] = D2[tt] + dD2
			if tt > 7/dt +1
				#I[tt+1] = I1[tt+1]+I2[tt+1]+R1[tt+1]+R2[tt+1]-(I1[tt-7/dt]+I2[tt-7/dt]+R1[tt-7/dt]+R2[tt-7/dt])
				#J1[tt+1] = (I1[tt+1]+R1[tt+1]+D1[tt+1]-D1[tt-7/dt]-(I1[tt-7/dt]+R1[tt-7/dt]))*100000/N1
				J2[tt+1] = (J1[tt+1] - J1[tt-Int(7/dt)])*100000/N1
				#J2[tt+1] = I2[tt+1]+R2[tt+1]+D2[tt+1]-D2[tt-7/dt]-(I2[tt-7/dt]+R2[tt-7/dt])
			end
		else
			print("boom")
			break
		end
	end
	return S1, I1, R1, D1, S2, I2, R2, D2, I, J1, J2
end;

# ╔═╡ 81f4c3f4-8e23-11eb-2f4e-f1b0173d3ef6
#Calculate the model
S1, I1, R1, D1, S2, I2, R2, D2, I, J1, J2 = sir_model(beta1,gamma1,delta1,rho1,N1,beta2,gamma2,delta2,rho2,N2,beta3,beta4,I0,J0,T,dt);

# ╔═╡ 116cc80c-8e25-11eb-1120-911f0eae6a3c
# Curves
plot(tt,[S1,I1,R1,D1],color=[:blue :red :green :yellow], label=["S" "I" "R" "D"], lw=5, grid=true,xlabel="Days",ylabel="Number of individuals")

# ╔═╡ c201ce30-8e26-11eb-15a1-99f683ff6d20
plot(tt,[S2,I2,R2,D2],color=[:blue :red :green :yellow], label=["S2" "I2" "R2" "D2"], lw=2, grid=true,xlabel="Days",ylabel="Number of individuals")

# ╔═╡ f7b43f62-8e27-11eb-1675-55e0b4f3462e
#plot incidence
plot(tt,J1,color=:blue,lw=2,xlabel="Days",ylabel="7 day incidence",label="I")

# ╔═╡ 5f5686ac-8e28-11eb-24a1-299e86f64808
#plot incidence
plot(tt,J2,color=:blue,lw=2,xlabel="Days",ylabel="7 day incidence",label="I")

# ╔═╡ Cell order:
# ╟─10a9374c-8e1a-11eb-10b4-fb8272dd46d1
# ╠═43970fa0-8e1b-11eb-3052-1701966a4478
# ╠═81f4c3f4-8e23-11eb-2f4e-f1b0173d3ef6
# ╠═883a4ad8-8e26-11eb-06c4-e30f3c0f70c1
# ╠═116cc80c-8e25-11eb-1120-911f0eae6a3c
# ╠═c201ce30-8e26-11eb-15a1-99f683ff6d20
# ╠═5f948160-8e27-11eb-1003-fbcf284800e6
# ╠═f7b43f62-8e27-11eb-1675-55e0b4f3462e
# ╠═5f5686ac-8e28-11eb-24a1-299e86f64808
# ╟─c30b3fb4-8e1c-11eb-05c7-d7f95aa97fcf
# ╠═01a236b0-8e1d-11eb-2630-ebff3f351bb3
# ╠═9f18ccd0-8e25-11eb-20e6-ab941f663550

max_stars_repo_name: roccminton/PlutoNotebooks

max_issues_repo_name: roccminton/PlutoNotebooks

max_forks_repo_name: roccminton/PlutoNotebooks


text:
con = createConn()

t = floor((now() - Dates.Month(6)), Dates.Day(1))

while t < Date("2018-03-09")

    #time now in my time zone

    startDate = floor((t - Dates.Day(14)), Dates.Minute(30))

    # Instrumen deklaration
    inst = instrument("RXH8 Comdty", "TRADE", startDate, t)

    #request data
    bundDaten = sort(BarFrame(con, inst, 30), cols = [:Volume])[(end-19):end,:]

    # Create File name
    path = "C:\\Users\\sduve\\JuliaPrPack\\bbgTS\\src\\"
    filename = string(path, "Bund", Date(t),".csv")

    # export File
    CSV.write(filename, bundDaten)

    t = DateTime(t) + Dates.Day(1)

end


text:
starting_numbers = open("starting_numbers.txt") do f
    parse.(Int, split(readline(f), ","))
end

function play_game(starting::Vector{Int64}, upto::Int64)::Int64
    last_seen = Dict{Int64, Tuple{Int64, Int64}}()
    prev = starting[end]

    for (i, v) = enumerate(starting)
        last_seen[v] = (i, 0)
        # println("Turn $(i): $(v)")
    end

    for i in (length(starting)+1):upto
        if last_seen[prev][2] == 0
            x = 0
        else
            x = i - 1 - last_seen[prev][2]
            # println("$(prev): $(i-1) - $(recent)")
        end
        # println("Turn $(i): $(x)")
        if haskey(last_seen, x)
            last_seen[x] = (i, last_seen[x][1])
        else
            last_seen[x] = (i, 0)
        end
        prev = x
    end

    return prev
end

# println(play_game([0,3,6], 10) == 0)
# println(play_game([0,3,6], 2020) == 436)
# println(play_game([1,3,2], 2020) == 1)
# println(play_game([2,1,3], 2020) == 10)
# println(play_game([1,2,3], 2020) == 27)
# println(play_game([2,3,1], 2020) == 78)
# println(play_game([3,2,1], 2020) == 438)
# println(play_game([3,1,2], 2020) == 1836)

@time println(play_game(starting_numbers, 2020))
@time println(play_game(starting_numbers, 30000000))


text:
"Constante universal dos gases em J / (kmol.K)"
const ℜ = 8314.472 #472 # Universal gas constant

"""
    myhorner(x, a, n)

Usa o método de Horner ara calcular o valor do polinômio
`y = a[1] + a[2]*x + a[3]*x^2 + ... + a[n]*x^(n-1)`

# Examples
```jldoctest
julia> a = [1, 2, 3];

julia> myhorner(0, a, 3)
1

julia> myhorner(1, a, 3)
6

julia> myhorner(2, a, 3)
17
```
"""
function myhorner(x, a, n)
    p = one(x) * a[n]

    for i = (n-1):-1:1
        p = muladd(x,p,a[i])
    end
    return p
end

"""
    IdealGas

Armazena um gás perfeito.
"""
struct IdealGas <: AbstractIdealGas
    "Massa molecular em kg/kmol"
    M::Float64
    "Coeficientes de polinômio representando o calor específico a pressão constante"
    cpcoefs::Vector{Float64}
end

"""
 Singleton to calculate Air properties.
"""
struct AirType <: AbstractIdealGas
end

"Massa molecular do Ar seco em kg/kmol"
const Mair = 28.9647 #kg/kmol

"Constante do gas para o ar"
const Rair = ℜ/Mair

"Constante representando o ar"
const Air = AirType()

"""
    heatcond(Air, 20.0)

Calcula a condutividade do ar em W/m.K

# Examples
```jldoctest
julia> heatcond(Air, 0.0)
0.02413198631696768

julia> heatcond(Air, 20.0)
0.025723538599749837

julia> heatcond(Air, 50.0)
0.028063733135454103
```
"""
function heatcond(Air::AirType, Tc=20.0)
    T = Tc + 273.15
    @evalpoly(T, 5.890e-02, 9.542e-02, -2.397e-05, -1.126e-08, 4.853e-12) * 1e-3
end

"""
    viscosity(Air, Tc)

Calcula a viscosidade dinâmica do ar em Pa.s usando a equação de Sutherland.

# Examples
```jldoctest
julia> viscosity(Air, 0.0)
1.7362296862138683e-5

julia> viscosity(Air, 20.0)
1.8369221075668256e-5

julia> viscosity(Air, 50.0)
1.982070679101198e-5
```
"""
function viscosity(Air::AirType, Tc)
    T = Tc + 273.15
    C = 120.0
    T₀ = 291.15
    μ₀ = 18.27e-6
    return μ₀ * (T₀ + C) / (T + C) * ( T/T₀ )^1.5
end

"""
    density(gas, Tc, Pk)

Calcula a densidade do ar em kg/m³ admitindo que o ar seja um gás perfeito. A temperatura `Tc` deve ser em °C e a pressão deve ser em kPa.

"""
function density(gas::IdealGas, Tc=20.0, Pk=101.325)
    return Pk*1000 * gas.M / (ℜ * (Tc + 273.15))
end

"""
    density(Air, Tc, Pk)

Calcula a densidade do ar em kg/m³ admitindo que o ar seja um gás perfeito. A temperatura `Tc` deve ser em °C e a pressão deve ser em kPa.

# Examples
```jldoctest
julia> density(Air, 0.0, 101.325)
1.2922595999236681

julia> density(Air, 20.0, 101.325)
1.2040958885183353

julia> density(Air, 20.0, 93)
1.105165730394327
```
"""
function density(Air::AirType, Tc=20.0, Pk=101.325)
    return Pk*1000 / (Rair * (Tc + 273.15))
end

"""
    kinviscosity(Air, Tc, Pk)

Calcula a viscosidade cinemática em m²/s do ar usando a equação de Sutherland.

# Examples
```jldoctest
julia> kinviscosity(Air, 20.0, 101.325)
1.5255613154091872e-5

julia> kinviscosity(Air, 0.0, 101.325)
1.3435610664578733e-5

julia> kinviscosity(Air, 20.0, 93)
1.6621236589659775e-5
```
"""
function kinviscosity(Air::AirType, Tc=20.0, P=101.325)
    return viscosity(Air, Tc) / density(Air, Tc, P)
end

"""
    specheat(Air, Tc)

Calcula o calor específico a pressão constante do ar em J/kg.K.

# Examples
```jldoctest
julia> specheat(Air, 0)
1003.6206378797417

julia> specheat(Air, 20)
1004.5273968591591
```
"""
function specheat(Air::AirType, Tc=20.0)
    T = Tc + 273.15
    return @evalpoly(T, 3.56839620E+00, -6.78729429E-04, 1.55371476E-06, -3.29937060E-12, -4.66395387E-13) * Rair  #, -1.06234659E+03,  3.71582965E+00) * Rair
end

"""
    thermaldiff(Air, Tc, Pk)

Calcula a difusividade térmica do ar seco em m²/s:

```math
α = k/(ρ.cₚ)
```

```jldoctest
julia> thermaldiff(Air, 0.0, 101.325)
1.8606886990841108e-5

julia> thermaldiff(Air, 20.0, 101.325)
2.1267079360835945e-5

julia> thermaldiff(Air, 20.0, 93)
2.3170825981039807e-5
```
"""
thermaldiff(Air::AirType, Tc=20.0, Pk=101.325) = heatcond(Air, Tc) / (density(Air, Tc, Pk) * specheat(Air, Tc))

"""
    prandtl(Air, Tc)

Número de Prandtl para o ar seco:

```math
Pr = ν/α = cₚμ/k
```

```jldoctest
julia> prandtl(Air, 0.0)
0.7220772970348109

julia> prandtl(Air, 20.0)
0.7173346605451432
```
"""
prandtl(Air::AirType, Tc=20.0) = viscosity(Air, Tc) * specheat(Air, Tc) / heatcond(Air, Tc)

"""
    volthermalexpansion(gas::AbstractIdealGas, Tc)

Coeficiente de expansão volumétrica em 1/K.
"""
volthermalexpansion(gas::AbstractIdealGas, Tc) = 1/(Tc + 273.15)


text:
# Portions translated from SLICOT-Reference distribution
# Copyright (c) 2002-2020 NICONET e.V.
function run_mc03md(datfile, io=stdout)
    NIN = 5
    NOUT = 6
    CP1MAX = 10
    CP2MAX = 10
    DP1MAX = 10
    DP2MAX = 10
    DP3MAX = 20
    RP1MAX = 10
    LDP11 = RP1MAX
    LDP12 = CP1MAX
    LDP21 = CP1MAX
    LDP22 = CP2MAX
    LDP31 = RP1MAX
    LDP32 = CP2MAX
    P1 = Array{Float64,3}(undef, LDP11,LDP12,DP1MAX+1)
    P2 = Array{Float64,3}(undef, LDP21,LDP22,DP2MAX+1)
    P3 = Array{Float64,3}(undef, LDP31,LDP32,DP3MAX+1)
    DWORK = Array{Float64,1}(undef, CP1MAX)
f = open(datfile,"r")
    readline(f)
    vs = split(readline(f))
    RP1 = parse(BlasInt, vs[1])
    CP1 = parse(BlasInt, vs[2])
    CP2 = parse(BlasInt, vs[3])
if ( RP1<0 || RP1>RP1MAX )
elseif ( CP1<0 || CP1>CP1MAX )
elseif ( CP2<0 || CP2>CP2MAX )
else
    vs = split(readline(f))
    DP1 = parse(BlasInt, vs[1])
if ( DP1<=-2 || DP1>DP1MAX )
else

    _isz,_jsz,_ksz = (RP1,CP1,DP1 + 1)
    for k in 1:_ksz
      for j in 1:_jsz
        vs = String[]
        while length(vs) < _isz
            append!(vs, replace.(split(readline(f)),'D'=>'E'))
        end
        P1[1:_isz,j,k] .= parsex.(Float64, vs)
      end
    end
    vs = split(readline(f))
    DP2 = parse(BlasInt, vs[1])
if ( DP2<=-2 || DP2>DP2MAX )
else

    _isz,_jsz,_ksz = (CP1,CP2,DP2 + 1)
    for k in 1:_ksz
      for j in 1:_jsz
        vs = String[]
        while length(vs) < _isz
            append!(vs, replace.(split(readline(f)),'D'=>'E'))
        end
        P2[1:_isz,j,k] .= parsex.(Float64, vs)
      end
    end
    vs = split(readline(f))
    DP3 = parse(BlasInt, vs[1])
if ( DP3<=-2 || DP3>DP3MAX )
else

    _isz,_jsz,_ksz = (RP1,CP2,DP3 + 1)
    for k in 1:_ksz
      for j in 1:_jsz
        vs = String[]
        while length(vs) < _isz
            append!(vs, replace.(split(readline(f)),'D'=>'E'))
        end
        P3[1:_isz,j,k] .= parsex.(Float64, vs)
      end
    end
    vs = split(readline(f))
    ALPHA = parse(Float64, replace(vs[1],'D'=>'E'))
# interp call 1

    DP3, INFO = SLICOT.mc03md!(RP1, CP1, CP2, DP1, DP2, DP3, ALPHA, P1, P2, P3)
    println(io, "DP3 = $DP3")
    @test INFO == 0
    INFO == 0 || return

if ( INFO!=0 )
else
    for i in 1:RP1
        for j in 1:CP2
            println(io, "elt($i,$j): ",P3[i,j,1:DP3+1])
        end
    end
end # if
end # if
end # if
end # if
end # if
    close(f)
end # run_mc03md()


text:
using ScottishTaxBenefitModel
using .FRSHouseholdGetter: initialise, get_household, get_num_households
using .STBParameters:
    TaxBenefitSystem,
    NationalInsuranceSys,
    IncomeTaxSys,
    weeklyise!,
    make_ubi_pre_adjustments!
import .Results: init_benefit_unit_result, BenefitUnitResult
using .ModelHousehold:
   BenefitUnit,
   Household,
   Person,
   child_pids,
   get_head,
   is_child,
   num_adults,
   num_children,
   make_eq_scales!

using .Definitions
import .ExampleHouseholdGetter

using .RunSettings: Settings

using .ExampleHelpers

using DataFrames,CSV, Dates

#
# full dataset is available ..
#
const IS_LOCAL = isdir("/mnt/data/frs/")

# pids for example people
# see make_pid
const RUK_PERSON = 320190010101

const SCOT_HEAD = 320190010201
const SCOT_SPOUSE = 320190010202

const DEFAULT_NUM_TYPE = Float64

const TEST_BASE_DATE = Date( 2021, 06, 14 ) # All the tests ran on this date; `TEST_BASE_DATE` is used
   # to adjust some ages e.g. social_security_ages_tests.jl so tests also pass next year.

function get_default_it_system(
   ;
  year     :: Integer=2019,
  scotland :: Bool = true,
  weekly   :: Bool = true )::Union{Nothing,IncomeTaxSys}
  it = nothing
  if year == 2019
     it = IncomeTaxSys{DEFAULT_NUM_TYPE}()
     if ! scotland
        it.non_savings_rates = [20.0,40.0,45.0]
        it.non_savings_thresholds = [37_500, 150_000.0]
        it.non_savings_basic_rate = 1
     end
     if weekly
        weeklyise!( it )
     end
  end
  it
end

function get_default_uc( ;
   year :: Integer = 2019,
   weekly :: Bool = true )
   uc = UniversalCreditSys{DEFAULT_NUM_TYPE}()
   if weekly
      weeklyise!( uc )
   else
      uc.taper /= 100
   end
   return uc
end

function to_nearest_p( x, y :: Real, ps :: Real = 1 ) :: Bool
   diff = abs(x-y)
   return diff <= 0.01*ps
   # round(x, digits=2) == round(y, digits=2)
end

function init_data(; reset :: Bool = false, settings = Settings() )
   nhh = get_num_households()
   num_people = -1
   if( nhh == 0 ) || reset
      @time nhh, num_people,nhh2 = initialise( settings )
   end
   (nhh,num_people)
end

function get_system(; scotland::Bool ) :: TaxBenefitSystem
   sys = TaxBenefitSystem{DEFAULT_NUM_TYPE}()
   weeklyise!(sys)
   # overwrite IT to get RuK system as needed
   # println( itn )
   sys.it = get_default_it_system( year=2019, scotland=scotland, weekly=true )
   #
   # So, assuming that the default has these things set ...
   #
   @assert sys.scottish_child_payment.amount > 0
   if ! scotland
      sys.scottish_child_payment.amount = 0.0
      sys.scottish_child_payment.maximum_age = 0
      sys.scottish_child_payment.qualifying_benefits = []
      sys.nmt_bens.carers.scottish_supplement = 0.0
   end
   if ! sys.ubi.abolished
      make_ubi_pre_adjustments!( sys )
   end
   return sys
end

function init_benefit_unit_result( bu :: BenefitUnit ) :: BenefitUnitResult
   return init_benefit_unit_result( DEFAULT_NUM_TYPE, bu )
end

#
# Messing around with periods to fit the Policy In Practice
# calculator.
#
const POLICY_IN_PRACTICE_WEEKS_PER_MONTH = 52/12
const PWPM = POLICY_IN_PRACTICE_WEEKS_PER_MONTH

"""
Compare PIP monthly calculations with our weekly, using their 52/12 weeks per month.
Throw an AssertionException if different, so we can use this directly
in a @test macro.
"""
function compare_w_2_m( uspw::Real, thempm::Real, ps :: Real = 1 ) :: Bool
   uspm = uspw*PWPM
   thempw = thempm/PWPM
   @assert to_nearest_p(uspw,thempw,ps) "us $(uspw)pw ($(uspm)pm) != $(thempm)pm"
   return true
end

"""
if the test was written at TEST_BASE_DATE, what age would we have to make somebody
to be sure the test will still work in some later year?
"""
function age_now( age :: Int ) :: Int
  yd = (Date(now()) - TEST_BASE_DATE).value ÷ 365 # leap years; no function for this
  return age + Int(yd)
end

max_stars_repo_name: grahamstark/ScottishTaxBenefitModel

max_issues_repo_name: grahamstark/ScottishTaxBenefitModel

max_forks_repo_name: grahamstark/ScottishTaxBenefitModel


text:
module TestSchema
type Post
    title::String
    author::String
end
end # module TestSchema

import Ecto
import Ecto.Schema: Assoc
using Base.Test

post = %(TestSchema.Post, author= "bar")
@test isa(post, Ecto.Schema.t)
@test Assoc([(:author, "bar")]) == post.struct


text:
using Unitful
using Printf

function shift_unit(u::U, d, n) where {U<:Unitful.Unit}
    i₀ = floor(Int, d/n)
    d = 3i₀
    iszero(d) && return u,0
    for (i,tens) in enumerate(Unitful.tens(u) .+ (d:(-3*sign(d)):0))
        haskey(Unitful.prefixdict, tens) && return U(tens, u.power),i-1 + i₀
    end
    u,0
end

power_step(::Any) = 3

function shift_unit(u::Unitful.FreeUnits, d, idx)
    tu = typeof(u)
    us,ds,a = tu.parameters

    uu,i = shift_unit(us[idx], d, power_step(1u))

    Unitful.FreeUnits{(uu,us[setdiff(1:end,idx)]...), ds, a}(),i
end

function si_round(q::Quantity, idx=1)
    v,u = ustrip(q), unit(q)
    if !iszero(v)
        u,i = shift_unit(u, log10(abs(v)), idx)
        q = u(q/(one(q)^i))
    end
     @sprintf("%.3g %s", ustrip(q), unit(q))
end


text:
@testset "valueselection.jl" begin

    include("../../RL/RL.jl")
    include("searchmetrics.jl")
    include("classic/basicheuristic.jl")
    include("learning/learning.jl")
end


text:
export DisjointSet
export find

mutable struct DisjointSet
    parent::Vector{Int64}

    DisjointSet(n::Int64) = new(-ones(Int64, n))
end

parent(set::DisjointSet, element::Int64) = set.parent[element]

function find(set::DisjointSet, element::Int64)
    p = parent(set, element)
    if p == -1 return element end
    find(set, p)
end

function Base.union!(set::DisjointSet, a::Int64, b::Int64)
    root1 = find(set, a)
    root2 = find(set, b)
    set.parent[root2] = root1
end


text:
using ZipFile

mutable struct NLModel <: NEOSModel
    solver::NEOSSolver
	xmlmodel::String
	last_results::String
    inner::AmplNLWriter.AmplNLMathProgModel
end
LinearQuadraticModel{S}(s::NEOSSolver{S, :NL}) = NLModel(
    s,
    "",
    "",
    AmplNLWriter.AmplNLMathProgModel("", String[], "")
)

function loadproblem!(m::NLModel, A::AbstractMatrix, x_l, x_u, c, g_l, g_u, sense)
    loadproblem!(AmplNLWriter.AmplNLLinearQuadraticModel(m.inner), A, x_l, x_u, c, g_l, g_u, sense)
end

function neos_writexmlmodel!(m::NLModel)
    # There is no non-linear binary type, only non-linear discrete, so make
    # sure binary vars have bounds in [0, 1]
    for i in 1:m.inner.nvar
        if m.inner.vartypes[i] == :Bin
            if m.inner.x_l[i] < 0
                m.inner.x_l[i] = 0
            end
            if m.inner.x_u[i] > 1
                m.inner.x_u[i] = 1
            end
        end
    end
    io = IOBuffer()
    print(io, "<model>")
    AmplNLWriter.make_var_index!(m.inner)
    AmplNLWriter.make_con_index!(m.inner)
    AmplNLWriter.write_nl_file(io, m.inner)
    print(io, "</model>")
    # Convert the model to NL and add
    m.xmlmodel = replace(m.solver.template, r"<model>.*</model>"is, String(take!(io)))
end

# Wrapper functions
for f in [:getvartype,:getsense,:status,:getsolution,:getobjval,:numvar,:numconstr,:get_solve_result,:get_solve_result_num,:get_solve_message,:get_solve_exitcode,:getsolvetime]
  @eval $f(m::NLModel) = $f(m.inner)
end
for f in [:setvartype!,:setsense!,:setwarmstart!]
  @eval $f(m::NLModel, x) = $f(m.inner, x)
end

function solfilename(job)
    "https://neos-server.org/neos/jobs/$(10_000 * floor(Int, job.number / 10_000))/$(job.number)-$(job.password)-ampl.sol"
end

getreducedcosts(m::NLModel) = fill(NaN, numvar(m))
getconstrduals(m::NLModel) = fill(NaN, numconstr(m))
getobjbound(m::NLModel) = NaN

function parseresults!(m::NLModel, job)
    if m.solver.print_level >= 2
        info("Getting solution file from $(solfilename(job))...")
    end
    res = get(solfilename(job))
	if res.status != 200
        error("Error retrieving results for job $(job.number):$(job.password). Response status is $(res.status).")
    end
    if m.solver.print_level >= 2
        info("Extracting file from .zip")
    end
    io = IOBuffer()
    write(io, res.data)
    if m.solver.print_level >= 2
        info("Reading results")
    end
    seekstart(io)
    AmplNLWriter.read_results(io, m.inner)
end


text:
"""
NonparametricRegression.jl implements local regression methods, also known as kernel regression.

Type `? npregress` for more detail.
"""
module NonparametricRegression

using LinearAlgebra
using Statistics
using StaticArrays
using DocStringExtensions

export NormalKernel, GaussianKernel, UniformKernel, EpanechnikovKernel, bandwidth
export localconstant, localconstantweights
export locallinear, llalphabeta, locallinearweights
export optimalbandwidth, leaveoneoutCV, optimizeAICc
export npregress

include("kernels.jl")
include("univariateopt.jl")

########################################
########################################
##### UTILITY FUNCTIONS
########################################
########################################

AICc(sigmasq,trH,Ny) = log(sigmasq) + 1 + 2*(trH + 1)/(Ny-trH-2)

## KDE
# this function is only used internally for trimmed LOOCV to prevent adding a dep on KernelDensity.jl
silvermanbw(x,alpha=0.9) = alpha * min(std(x), only(diff(quantile(x, [0.25, 0.75])))/1.34) * length(x)^(-1/5)

function densityestimator(x,xgrid=x; h=silvermanbw(x), kernelfun=NormalKernel)
	K = kernelfun(h)
	N = length(x)
	L = length(xgrid)
	f = zeros(L) 	# preallocate
	@inbounds for i in 1:L
		for j in 1:N
			f[i] += unscaledkern(K,x[j],xgrid[i])
		end
		f[i] *= scaling(K) / N
	end
	return f
end

########################################
########################################
##### SPECIFIC NONPARAMETRIC REGRESSION METHODS
########################################
########################################

## LOCAL CONSTANT aka NADARAYA-WATSON ESTIMATOR
function localconstant(x,y,xgrid, h; kernelfun=NormalKernel)
	K = kernelfun(h)
	Nx = length(x)
	@assert Nx == length(y) "Length of `x` and `y` must agree."
	L = length(xgrid)
	m = zeros(L)	# m is expected y, at points xgrid
	@inbounds for i in 1:L
		wu = 0.0
		wsum = 0.0
		for j in 1:Nx
			ku = unscaledkern(K,x[j],xgrid[i])
			wu += ku * y[j]
			wsum += ku
		end
		m[i] = wu / wsum
	end
	return m
end

function localconstant_aicc(x,y,h; kernelfun=NormalKernel)
	K = kernelfun(h)
	Nx = length(x)
	@assert Nx == length(y) "Length of `x` and `y` must agree."

	sigmasq = 0.0
	traceH = 0.0
	@inbounds for i in 1:Nx
		wu = 0.0
		wsum = 0.0
		tH = 0.0
		for j in 1:Nx
			ku = unscaledkern(K,x[j],x[i])
			wu += ku * y[j]
			wsum += ku
			if j == i
				tH = ku
			end
		end
		yp = wu / wsum
		sigmasq += (y[i] - yp)^2 / Nx
		traceH += tH / wsum
	end
	aicc = AICc(sigmasq,traceH,Nx)
	return aicc
end

function localconstantweights(x,xgrid, h; kernelfun=NormalKernel, normalizeweights::Bool=true)
	K = kernelfun(h)
	Nx = length(x)
	L = length(xgrid)
	W = zeros(Nx,L)	# m is expected y, at points xgrid
	@inbounds for i in 1:L
		wsum = 0.0
		for j in 1:Nx
			ku = unscaledkern(K,x[j],xgrid[i])
			W[j,i] = ku
			wsum += ku
		end
		if normalizeweights
			W[:,i] ./= wsum
		end
	end
	return W
end

## LOCAL LINEAR REGRESSION
function locallinear(x,y,xgrid, h; kernelfun=NormalKernel)
	K = kernelfun(h)
	Nx = length(x)
	@assert Nx == length(y) "Length of `x` and `y` must agree."
	L = length(xgrid)
	m = zeros(L)	# m is expected y, at points xgrid
	@inbounds for i in 1:L
		s1sum = 0.0
		s2sum = 0.0
		wsum = 0.0
		for j in 1:Nx
			ku = unscaledkern(K,x[j],xgrid[i])
			xmx = x[j]- xgrid[i]
			s1sum += ku*xmx
			s2sum += ku*xmx^2
			wsum += ku
		end
		s1sum /= wsum
		s2sum /= wsum

		wu = 0.0
		wsum = 0.0
		for j in 1:Nx
			ku = unscaledkern(K,x[j],xgrid[i])
			xmx = x[j]- xgrid[i]
			w = ku * (s2sum - s1sum*xmx)
			wu += w * y[j]
			wsum += w
		end
		m[i] = wu / wsum
	end
	return m
end

function locallinear_aicc(x,y, h; kernelfun=NormalKernel)
	K = kernelfun(h)
	Nx = length(x)
	@assert Nx == length(y) "Length of `x` and `y` must agree."

	sigmasq = 0.0
	traceH = 0.0
	@inbounds for i in 1:Nx
		s1sum = 0.0
		s2sum = 0.0
		wsum = 0.0
		for j in 1:Nx
			ku = unscaledkern(K,x[j],x[i])
			xmx = x[j]- x[i]
			s1sum += ku*xmx
			s2sum += ku*xmx^2
			wsum += ku
		end
		s1sum /= wsum
		s2sum /= wsum

		wu = 0.0
		wsum = 0.0
		tH = 0.0
		for j in 1:Nx
			ku = unscaledkern(K,x[j],x[i])
			xmx = x[j]- x[i]
			w = ku * (s2sum - s1sum*xmx)
			wu += w * y[j]
			wsum += w
			if j == i
				tH = w
			end
		end
		yp = wu / wsum
		sigmasq += (y[i] - yp)^2 / Nx
		traceH += tH / wsum
	end
	aicc = AICc(sigmasq,traceH,Nx)
	return aicc
end

function llalphabeta(x,y,xgrid, h; kernelfun=NormalKernel)
	K = kernelfun(h)
	Nx = length(x)
	@assert Nx == length(y) "Length of `x` and `y` must agree."
	L = length(xgrid)
	avec = zeros(L)
	bvec = zeros(L)
	@inbounds for i in 1:L
		l = SMatrix{2,2}(0.0,0.0,0.0,0.0)
		r = SVector(0.0,0.0)
		wsum = 0.0
		for j in 1:Nx
			ku = unscaledkern(K,x[j],xgrid[i])
			zj = SVector(1.0,x[j]- xgrid[i])
			l += ku*zj*zj'
			r += ku*zj*y[j]
			wsum += ku
		end
		a,b = l/wsum \ r/wsum
		avec[i] = a
		bvec[i] = b
	end
	return avec,bvec
end

function locallinearweights(x,xgrid, h; kernelfun=NormalKernel, normalizeweights::Bool=true)
	K = kernelfun(h)
	Nx = length(x)
	L = length(xgrid)
	W = zeros(Nx,L)
	for i in 1:L
		s1sum = 0.0
		s2sum = 0.0
		wsum = 0.0
		for j in 1:Nx
			ku = unscaledkern(K,x[j],xgrid[i])
			xmx = x[j]- xgrid[i]
			s1sum += ku*xmx
			s2sum += ku*xmx^2
			wsum += ku
		end
		s1sum /= wsum
		s2sum /= wsum

		wsum = 0.0
		for j in 1:Nx
			ku = unscaledkern(K,x[j],xgrid[i])
			xmx = x[j]- xgrid[i]
			w = ku * (s2sum - s1sum*xmx)
			W[j,i] = w
			wsum += w
		end
		if normalizeweights
			W[:,i] ./= wsum
		end
	end
	return W
end

########################################
########################################
##### BANDWIDTH SELECTION
########################################
########################################
## LOOCV
function leaveoneoutCV_mse(h,x,y; kernelfun=NormalKernel, method=:lc, trimmed=true)
	if trimmed
		kde = densityestimator(x)
	end

	Nx = length(x)
	@assert Nx == length(y) "Length of `x` and `y` must agree."
	mse = 0.0
	inds = eachindex(x)
	@inbounds for i in 1:Nx
		xtmp = view(x,filter(!=(i),inds))
		ytmp = view(y,filter(!=(i),inds))
		# pred = estimatorfun(xtmp,ytmp,x[i], h, kernelfun)
		if method in (:nw,:lc,:localconstant)
			pred = only(localconstant(xtmp,ytmp,x[i], h; kernelfun))
		elseif method in (:ll,:locallinear,:llalphabeta)
			pred = only(first(llalphabeta(xtmp,ytmp,x[i], h; kernelfun)))
		else
			error("Unknown method supplied")
		end
		if trimmed
			mse += (y[i] - pred)^2 * kde[i] / Nx
		else
			mse += (y[i] - pred)^2 / Nx
		end
	end

	return mse
end

function leaveoneoutCV(x, y; kernelfun=NormalKernel, method=:lc, hLB = silvermanbw(y)/100, hUB = silvermanbw(y)*100,trimmed=true)
	objfun(h) = leaveoneoutCV_mse(h,x,y; kernelfun, method, trimmed)
	opt = optimize(objfun,hLB,hUB)
	# @assert opt.converged "Convergence failed, cannot find optimal bandwidth."
	return opt
end

function estimatorAICc(h, x, y; kernelfun=NormalKernel, method=:lc)
    if method in (:nw,:lc,:localconstant)
		return localconstant_aicc(x, y, h; kernelfun)
	elseif method in (:ll,:locallinear,:llalphabeta)
		return locallinear_aicc(x, y, h; kernelfun)
	else
		error("Unknown method supplied")
	end
end

function optimizeAICc(x, y; kernelfun=NormalKernel, method=:lc, hLB = silvermanbw(y)/100, hUB = silvermanbw(y)*100)
    objfun(h) = estimatorAICc(h,x,y; kernelfun, method)
	opt = optimize(objfun,hLB,hUB)
	# @assert opt.converged "Convergence failed, cannot find optimal bandwidth."
	return opt
end

########################################
########################################
##### CONVENIENCE METHODS
########################################
########################################
"""
$SIGNATURES

Search for the optimal bandwidth to use for the local regression of `y` against `x` using `method ∈ (:lc,:ll)`.

The keyword argument `bandwidthselection` should be `:aicc` for the bias-corrected AICc method or `:loocv` for leave-one-out cross validation. `hLB` and `hUB` are the lower and upper bounds used when searching for the optimal bandwidth.
"""
function optimalbandwidth(x, y; kernelfun=NormalKernel, method=:lc, bandwidthselection=:aicc, hLB = silvermanbw(y)/100, hUB = silvermanbw(y)*100)
    # automatically select bandwidth
    if bandwidthselection in (:cv, :loocv, :leaveoneoutcv)
        h = leaveoneoutCV(x,y; kernelfun, method, hLB, hUB)
    elseif bandwidthselection in (:aicc, :aic)
        h = optimizeAICc(x,y; kernelfun, method, hLB, hUB)
    else
        error("Unknown bandwidth selection procedure supplied. See `? optimalbandwidth` for valid inputs.")
    end
    return h
end

"""
$SIGNATURES

Estimate a local regression of `y` against `x` evaluated at the values `xgrid` using bandwidth `h`.

The keyword argument `method` can be `:lc` for a local constant estimator (Nadaraya-Watson) or `:ll` for a local linear estimator.

The keyword argument `kernelfun` should be a function that constructs a kernel with a given bandwidth. Defaults to `NormalKernel` (defined and exported by this package.)
"""
function npregress(x,y,xgrid,h; kernelfun=NormalKernel, method=:lc)
    if method in (:nw,:lc,:localconstant)
		return localconstant(x,y,xgrid,h; kernelfun)
	elseif method in (:ll,:locallinear)
		return locallinear(x,y,xgrid,h; kernelfun)
	else
		error("Unknown regression method supplied. See `? npregress` for valid inputs.")
	end
end

"""
$SIGNATURES

Estimate a local regression of `y` against `x`, possibly evaluated at the values `xgrid`, and automatically select the optimal bandwidth.

The keyword argument `method` can be `:lc` for a local constant estimator (Nadaraya-Watson) or `:ll` for a local linear estimator.

The keyword argument `bandwidthselection` should be `:aicc` for the bias-correct AICc method or `:loocv` for leave-one-out cross validation. `hLB` and `hUB` are the lower and upper bounds used when searching for the optimal bandwidth.

The keyword argument `kernelfun` should be a function that constructs a kernel with a given bandwidth. Defaults to `NormalKernel` (defined and exported by this package.)
"""
function npregress(x, y, xgrid=x; kernelfun=NormalKernel, method=:lc, bandwidthselection=:aicc, hLB = silvermanbw(y)/100, hUB = silvermanbw(y)*100)
    h = optimalbandwidth(x,y; kernelfun, method, bandwidthselection, hLB, hUB)
    return npregress(x,y,xgrid,h; kernelfun, method)
end

########################################
########################################
##### CURRENTLY UNUSED METHODS (BUT MAYBE RELEVANT FOR FUTURE WORK)
########################################
########################################
## commented out to increase coverage :)

# function transformedresiduals(y::AbstractVector{T1},H::AbstractMatrix{T2}; method=:HC3) where {T1,T2}
# 	r = y .- H*y
# 	if method == :HC0
# 		return r
# 	elseif method == :HC1
# 		error("HC1 not implemented")
# 	elseif method == :HC2
# 		Hd = diag(H)
# 		return r./sqrt.(1 .- Hd)
# 	elseif method == :HC3
# 		Hd = diag(H)
# 		return r./(1 .- Hd)
# 	else
# 		error("unknown method")
# 	end
# end

end # module


text:
@testset "parse single subject" begin
    bids_root = joinpath(
        @__DIR__,
        "data",
        "bids_root"
    )

    my_layout = Layout(bids_root)

    @test my_layout.root == bids_root
    @test my_layout.longitudinal == true
    @test my_layout.description == OrderedDict{String,Any}()
    @test my_layout.subjects_detail == DataFrame()
    @test total_subjects(my_layout) == 1
    @test total_sessions(my_layout) == 1
    @test total_files(my_layout) == 1
    @test sprint(show, my_layout) == """
        Layout:
            root = $bids_root
            total subject = 1
            total session = 1
            total files = 1"""
end

@testset "query files" begin
    # get_files in layout implements all of the get_files (from file subject, and
    # session), so testing the others are not necessary
    bids_root = joinpath(
        @__DIR__,
        "data",
        "bids_root"
    )

    my_layout = Layout(bids_root)

    f = get_files(my_layout, key1="value1", run="001")
    @test length(f) == 1
    @test f[1].path == joinpath(
        @__DIR__,
        "data",
        "bids_root",
        "sub-subtest",
        "ses-1",
        "test",
        "sub-test_ses-1_run-001_modlbl.nii.gz"
    )

    f = get_files(my_layout, key1="value1", run="002")
    @test length(f) == 0

    f = get_files(my_layout, key1="value2", run="001")
    @test length(f) == 0

    f = get_files(my_layout, random_key="value1")
    @test length(f) == 0

    # Test with path argument
    f = get_files(my_layout, path=".nii.gz", key1="value1", run="001")
    @test length(f) == 1

    f = get_files(my_layout, path=".json", key1="value1", run="001")
    @test length(f) == 0
end


text:
#**************************************************************************************
# Physics_Model.jl
# ================ part of the GeoEfficiency.jl package.
#
# here is the place where all physical elements is being modeled and created as computer objects.
#
#**************************************************************************************

#------------------ consts - globals - imports -------------------

import Base: show, isless

#-------------------------- Point ----------------------------------

"""

    Point(Height::Real, Rho::Real)

construct and return a `Point` source.
The `Point` can be used as either a source by itself or an `anchor point` of a higher dimension source.

*  `Height` : point height relative to the detector surface.
*  `Rho` : point off-axis relative to the detector axis of symmetry.

!!! warning "Interpretation of `Height`"
    Each detector type give different interpretation to the `Height` as follow:-
    *  for `CylDetector` the point source `height` is consider to be measured
       from the detector `face surface`.
    *  for `BoreDetector` the point source `height` is consider to be measured
       from the `detector middle`, +ve value are above the detector center while -ve are below.
    *  for `WellDetector` the point source `height` is considered to be measured
       from the detector `hole surface`.

"""
struct Point
   	Height::Float64
   	Rho::Float64
   	Point(Height::Float64, Rho::Float64) = new(Height, Rho)
end #type
Point(Height::Real, Rho::Real) = Point(float(Height), float(Rho))

"""

	Point(Height::Real)

construct and return an `axial point`.

**see also:** [`Point(Height::Real, Rho::Real)`](@ref).

"""
Point(Height::Real) = Point(Height, 0.0)

"""

	Point()

construct and return a `point`. prompt to input information via the `console`.

**see also:** [`Point(Height::Real, Rho::Real)`](@ref).

"""
function Point()
   	printstyled("\n II- The Radioactive Source Anchor Point:-\n", color = :yellow)
   	Height = getfloat("\n\t > Height (cm) = ")
   	Rho = getfloat("\n\t > Off-axis (cm) = ")
   	Point(Height, Rho)
end #function

"""

	Point(xHeight::Real, aPnt::Point)

construct and return a `point` that has the same off-axis distance as `aPnt` but of new
height `xHeight`.

**see also:** [`Point(Height::Real, Rho::Real)`](@ref)

"""
Point(xHeight::Real, aPnt::Point) = Point(xHeight, aPnt.Rho)

"""

	Point(aPnt::Point, xRho::Real)

construct and return a `point` that has the same height as `aPnt` but of new
off-axis distance `Rho`.

**see also:** [`Point(Height::Real, Rho::Real)`](@ref).

"""
Point(aPnt::Point, xRho::Real) = Point(aPnt.Height, xRho)

id(aPnt::Point) = "Point[Height=$(aPnt.Height), Rho=$(aPnt.Rho)]"
show(pnt::Point) = print(id(pnt))

#--------------source---------------------------------------------

"""

	source(anchorPnt::Point = Point())

return a tuple that describe the source (`anchorPnt`, `SrcRadius`, `SrcLength`) according to
the input from the `console`.

*  `anchorPnt` : the source anchoring point. if it is missing the user is prompt
   to input it via the `console`.
*  `SrcRadius` : source radius.
*  `SrcLength` : source length.

!!! warning "Point/Cylinder Source"
	if source type set to point source, both `SrcRadius` and `SrcLength` are set to zero.
    for more information **see also:** [`typeofSrc()`](@ref) and [`typeofSrc(x::Int)`](@ref).
"""
function source(anchorPnt::Point = Point())::Tuple{Point,Float64,Float64}

   	if setSrcToPoint()
      		@warn("""srcType is set to ``srcPoint``,
		**see** `setSrcToPoint` for more information.""", _file=nothing)
      		return (anchorPnt, 0.0, 0.0)
   	end #if

   	SrcRadius = getfloat("\n\t > Source Radius (cm) = ", 0.0)
    if 0.0 != SrcRadius
        SrcLength = getfloat("\n\t > Source Length (cm) = ", 0.0)
      	@error("currently only axial non-point sources are allowed", _file=nothing)
      	@warn("the off-axis will be set to `Zero`", _file=nothing)
      	anchorPnt = Point(anchorPnt, 0.0)

   	else
        SrcLength = 0.0
      	@warn("`SrcLength` is set to `zero`", _file=nothing)

   	end #if
    return (anchorPnt, SrcRadius, SrcLength)
end #function

#---------------- Detector ------------------------------------

"""
abstract super-supertype of all detectors types
"""
abstract type RadiationDetector end

"""

	Detector

abstract supertype of all detectors types of cylidericalish shapes.
also can be used to construct any leaf type.

"""
abstract type Detector <: RadiationDetector end
show(io::IO, detector::RadiationDetector) = print(id(detector))
isless(detector1::RadiationDetector, detector2::RadiationDetector) = isless(volume(detector1), volume(detector2))

##-------------- CylDetector -----------------------------------

"""

	CylDetector(CryRadius::Real, CryLength::Real)

construct and return a `cylindrical` detector of the given crystal dimensions:-

*  `CryRadius` : the detector crystal radius.
*  `CryLength` : the detector crystal length.

!!! warning "Invalid Arguments"
    both `CryRadius` and `CryLength` should be `positive`, while `CryLength` can also be set to **`zero`**.

"""
struct CylDetector <: Detector
   	CryRadius::Float64    	#Real
    CryLength::Float64		#Real

   	function CylDetector(CryRadius::Float64, CryLength::Float64)
      		@validateDetector	CryRadius > 0.0		"Crystal Radius: expect +ve number, get '$(CryRadius)'"
      		@validateDetector	CryLength >= 0.0	"Crystal Length: expect +ve number or zero, get '$(CryLength)'"
      		new(CryRadius, CryLength)
   	end #function

end #type
CylDetector(CryRadius::Real, CryLength::Real) = CylDetector(float(CryRadius), float(CryLength))

"""

    CylDetector(CryRadius::Real)

construct and return a `cylindrical` (really `disk`) detector with crystal length equal to **`zero`**.

**see also:** [`CylDetector(CryRadius::Real, CryLength::Real)`](@ref).

"""
CylDetector(CryRadius::Real) = CylDetector(CryRadius, 0.0)

"""

    CylDetector()

construct and return a `cylindrical` detector according to the input from the `console`.

**see also:** [`CylDetector(CryRadius::Real, CryLength::Real)`](@ref).

"""
function CylDetector()
   	printstyled(" I- The Cylindrical Detector physical Dimensions:-\n", color = :yellow)
   	CryRadius = getfloat("\n\t > Crystal Radius (cm) = ", 0.0, lower=false)
   	CryLength = getfloat("\n\t > Crystal Length (cm) = ", 0.0)
   	CylDetector(CryRadius, CryLength)
end #function

id(detector::CylDetector) = "CylDetector[CryRadius=$(detector.CryRadius), CryLength=$(detector.CryLength)]"
volume(detector::CylDetector) = pi * detector.CryRadius^2 * detector.CryLength

##------------- BoreDetector -------------------------------------

"""

	BoreDetector(CryRadius::Real, CryLength::Real, HoleRadius::Real)

construct and return a `bore-hole` detector of the given crystal dimensions:-

*  `CryRadius` : the detector crystal radius.
*  `CryLength` : the detector crystal length.
*  `HoleRadius` : the detector hole radius.

!!! warning "Invalid Arguments"
    `CryRadius` and `CryLength`, `HoleRadius` should be `positive` numbers, also
    `CryRadius` should be greater than `HoleRadius`.

"""
struct BoreDetector <: Detector
   	CryRadius::Float64    	#Real
    CryLength::Float64    	#Real
   	HoleRadius::Float64    	#Real

   	function BoreDetector(CryRadius::Float64, CryLength::Float64, HoleRadius::Float64)
      		@validateDetector	CryRadius > 0.0		"Crystal Radius: expect +ve number, get '$(CryRadius)'"
      		@validateDetector	CryLength > 0.0		"Crystal Length: expect +ve number, get '$(CryLength)'"
      		@validateDetector	CryRadius > HoleRadius > 0.0	"Hole Radius: expect +ve number Less than 'Crystal Radius=$(CryRadius)', get $(HoleRadius)."
      		new(CryRadius, CryLength, HoleRadius)
   	end #function

end #type
BoreDetector(CryRadius::Real, CryLength::Real, HoleRadius::Real) = BoreDetector(float(CryRadius), float(CryLength), float(HoleRadius))

"""

	BoreDetector()

construct and return a `bore-hole` detector according to the input from the `console`.

**see also:** [`BoreDetector(CryRadius::Real, CryLength::Real, HoleRadius::Real)`](@ref).

"""
function BoreDetector()
   	printstyled(" I- The Bore Hole Detector physical Dimensions:-\n", color = :yellow)
   	CryRadius  = getfloat("\n\t > Crystal Radius (cm) = ", 0.0, lower=false)
   	CryLength  = getfloat("\n\t > Crystal Length (cm) = ", 0.0, lower=false)
   	HoleRadius = getfloat("\n\t > Hole Radius (cm) = ", 0.0, CryRadius)
   	BoreDetector(CryRadius, CryLength, HoleRadius)
end #function

id(detector::BoreDetector) = "BoreDetector[CryRadius=$(detector.CryRadius), CryLength=$(detector.CryLength), HoleRadius=$(detector.HoleRadius)]"
volume(detector::BoreDetector) = pi * (detector.CryRadius^2 - detector.HoleRadius^2) * detector.CryLength

##----------------------- WellDetector ------------------------------------------

"""

	WellDetector(CryRadius::Real, CryLength::Real, HoleRadius::Real, HoleDepth::Real)

construct and return a `Well-Type` detector of the given crystal dimensions:-

*  `CryRadius` : the detector crystal radius.
*  `CryLength` : the detector crystal length.
*  `HoleRadius` : the detector hole radius.
*  `HoleDepth` : the detector hole length.

!!! warning "Invalid Arguments"
    all arguments should be `positive` numbers, also
    `CryRadius` should be greater than `HoleRadius` and
    `CryLength` should be greater than  `HoleDepth`.

"""
struct WellDetector <: Detector
   	CryRadius::Float64
    CryLength::Float64
   	HoleRadius::Float64
   	HoleDepth::Float64

   	function WellDetector(CryRadius::Float64, CryLength::Float64, HoleRadius::Float64, HoleDepth::Float64)
      		@validateDetector	CryRadius > 0.0				"Crystal Radius: expect +ve number, get '$(CryRadius)'"
      		@validateDetector	CryLength > 0.0				"Crystal Length: expect +ve number, get '$(CryLength)'"
      		@validateDetector	CryRadius > HoleRadius > 0.0	"Hole Radius: expect +ve number less than 'Crystal Radius = $(CryRadius)', get '$(HoleRadius)'"
      		@validateDetector	CryLength > HoleDepth > 0.0	   	"Hole Depth: expect +ve number less than 'Crystal Length = $(CryLength)', get '$(HoleDepth)'"
      		new(CryRadius, CryLength, HoleRadius, HoleDepth)
   	end #if

end #type
WellDetector(CryRadius::Real, CryLength::Real, HoleRadius::Real, HoleDepth::Real) = WellDetector(float(CryRadius), float(CryLength), float(HoleRadius), float(HoleDepth))

"""

	WellDetector()

construct and return a Well-Type detector according to the input from the `console`.

**see also:** [`WellDetector(CryRadius::Real, CryLength::Real, HoleRadius::Real, HoleDepth::Real)`](@ref).

"""
function WellDetector()
   	printstyled(" I- The Well-Type Detector physical Dimensions:-\n", color = :yellow)
   	CryRadius  = getfloat("\n\t > Crystal Radius (cm) = ", 0.0, lower=false)
   	CryLength  = getfloat("\n\t > Crystal Length (cm) = ", 0.0, lower=false)
   	HoleRadius = getfloat("\n\t > Hole Radius (cm) = ", 0.0, CryRadius)
   	HoleDepth  = getfloat("\n\t > Hole Radius (cm) = ", 0.0, CryLength, lower=false)
   	WellDetector(CryRadius, CryLength, HoleRadius, HoleDepth)
end #function

id(detector::WellDetector) = "WellDetector[CryRadius=$(detector.CryRadius), CryLength=$(detector.CryLength), HoleRadius=$(detector.HoleRadius), HoleDepth=$(detector.HoleDepth)]"
volume(detector::WellDetector) = pi * (detector.CryRadius^2 * detector.CryLength - detector.HoleRadius^2 * detector.HoleDepth)

#---------------------------- Detector -------------------------------------

"""

	Detector()

construct and return an object of the `Detector` leaf types
(`CylDetector`, `BoreDetector` or `WellDetector`) according to the input from the console.

!!! note
    all required information is acquired from the `console` and would warn user on invalid data.

"""
function Detector()
   	printstyled("\n I- The detector physical Dimensions :-\n", color = :yellow)
   	CryRadius  = getfloat("\n\t > Crystal Radius (cm) = ", 0.0, lower=false)
   	CryLength  = getfloat("\n\t > Crystal Length (cm) = ", 0.0)
   	HoleRadius = getfloat("\n(zero for cylindrical detectors) > Hole Radius (cm) = ", 0.0, CryRadius, upper=false)
   	if   0.0 == HoleRadius
      		return CylDetector(CryRadius, CryLength)

   	else
      		HoleDepth = getfloat("\n(zero for Bore-Hole detectors) > Hole Depth (cm) = ", 0.0, CryLength)
      		return 0.0 == HoleDepth ?
				BoreDetector(CryRadius, CryLength, HoleRadius) :
				WellDetector(CryRadius, CryLength, HoleRadius, HoleDepth)
   	end #if
end #function

"""

	Detector(CryRadius::Real)

same as [`CylDetector(CryRadius::Real)`](@ref).

"""
Detector(CryRadius::Real) = CylDetector(CryRadius)

"""

	Detector(CryRadius::Real, CryLength::Real)

same as [`CylDetector(CryRadius::Real, CryLength::Real)`](@ref).

"""
Detector(CryRadius::Real, CryLength::Real) = CylDetector(CryRadius, CryLength)

"""

	Detector(CryRadius::Real, CryLength::Real, HoleRadius::Real)

same as [`BoreDetector(CryRadius::Real, CryLength::Real, HoleRadius::Real)`](@ref) except
when `HoleRadius` = `0.0` it acts as  [`CylDetector(CryRadius::Real, CryLength::Real)`](@ref).

"""
Detector(CryRadius::Real, CryLength::Real, HoleRadius::Real) = 0.0 == HoleRadius ?
				                                        CylDetector(CryRadius, CryLength) :
				                                        BoreDetector(CryRadius, CryLength, HoleRadius)

"""

	Detector(CryRadius::Real, CryLength::Real, HoleRadius::Real, HoleDepth::Real)

construct and return `well-type`, `bore-hole` or `cylindrical` detector according to the arguments.
it inspect the arguments and call the appropriate leaf type constructor.

!!! note "Missing Argument(s)"
    if the value(s) of the last argument(s) is\are `zero`, it acts as a missing argument(s).

**see also:** [`CylDetector`](@ref), [`BoreDetector`](@ref), [`WellDetector`](@ref).

"""
Detector(CryRadius::Real, CryLength::Real, HoleRadius::Real, HoleDepth::Real) = 0.0 == HoleDepth ?
				                                        Detector(CryRadius, CryLength, HoleRadius) :
				                                        WellDetector(CryRadius, CryLength, HoleRadius, HoleDepth)

"""

	Detector(detector::Detector)

return the inputted detector.

"""
Detector(detector::Detector) = detector

"""

	Detector(detectors::Vector{<: Detector})

convert the array `detectors` of any of the leaf `Detector` types
to an array of type `Detector`.

"""
Detector(detectors::Vector{<:Detector}) = Detector[detectors...]


text:
# License for this file: MIT (expat)
# Copyright 2017-2018, DLR Institute of System Dynamics and Control
#
# This file is part of module
#   Modia3D.DLR_Visualization (Modia3D/renderer/DLR_Visualization/_module.jl)
#

getVisualElement(data)               = data
getVisualElement(data::Solids.Solid) = data.geo

# Get function to visualize AbstractObject3Ddata object:
# - element  = getVisualElement(data): Get element from AbstractObject3Ddata that shall be visualized
# - name     = Basics.trailingPartOfTypeAsString( element ): Get trailing solid of typeof(element) as String (e.g. "Sphere")
# - symbol   = Symbol("visualize", name): Construct function name as symbol (e.g. :visualizeSphere)
# - function = getfield(DLR_Visualization, symbol): Get a Function pointer to symbol stored in module DLR_Visualization
#
# An alternative and even more efficient method would be to determine the C-function pointer with "cfunction(..)"
# and then call this function with "ccall(..)". However, if something goes wrong, Julia will crash.
# Therefore, this solution is not (yet) used.
getVisualizeFunction(data) = getfield(DLR_Visualization, Symbol("visualize", Basics.trailingPartOfTypeAsString( getVisualElement(data) )))

function Composition.initializeVisualization(renderer::Modia3D.AbstractDLR_VisualizationRenderer, velements::Vector{Composition.Object3D})::NOTHING
  simVis = renderer.simVis
  @assert(length(velements) > 0)
  if simVis.isInitialized
     Composition.closeVisualization(renderer)
  end

  # Initialize SimVis
  SimVis_init(simVis)

  # Determine visualize functions and ids of visual elements
  for i in eachindex(velements)
    push!(simVis.visualize, getVisualizeFunction(velements[i].data))
    push!(simVis.ids      , SimVis_getObjectID(simVis,0))
  end
  simVis.velements = velements
  simVis.isInitialized = true
  return nothing
end

function Composition.visualize!(renderer::Modia3D.AbstractDLR_VisualizationRenderer, time::Float64)
  simVis = renderer.simVis
  @assert(length(simVis.velements) > 0)
  if simVis.isInitialized
     velements = simVis.velements
     visualize = simVis.visualize
     ids       = simVis.ids
     for i in eachindex(velements)
        visualize[i](velements[i].data, velements[i], ids[i], simVis)
     end
     SimVis_setTime(simVis, time)
  else
     error("visualize! called without initializing visualization first.")
  end
end

function Composition.closeVisualization(renderer::Modia3D.AbstractDLR_VisualizationRenderer)
  simVis = renderer.simVis
  if simVis.isInitialized
    sleep(0.2)
    for id in simVis.ids
      SimVis_freeObjectID(simVis, id)
    end
    SimVis_shutdown(simVis)
  end
  empty!(simVis.ids)
  empty!(simVis.visualize)
  empty!(simVis.velements)
  simVis.isInitialized = false
end

Composition.isVisible(data::Modia3D.AbstractVisualElement, renderer::ProfessionalEdition) = true
Composition.isVisible(data::Solids.Solid                 , renderer::ProfessionalEdition) = typeof(data.material) != NOTHING && typeof(data.geo) != NOTHING

Composition.isVisible(data::Graphics.Spring              , renderer::CommunityEdition) = true
Composition.isVisible(data::Graphics.GearWheel           , renderer::CommunityEdition) = true
Composition.isVisible(data::Graphics.CoordinateSystem    , renderer::CommunityEdition) = true
Composition.isVisible(data::Modia3D.AbstractGeometry     , renderer::CommunityEdition) = true
Composition.isVisible(data::Solids.Solid                 , renderer::CommunityEdition) = typeof(data.material) != NOTHING && typeof(data.geo) != NOTHING


text:
struct VirtualGate{T<:AbstractFloat} <: Wall{T}
    sp::SVector{2,T}
    ep::SVector{2,T}
    normal::SVector{2,T}
    width::T
    center::SVector{2,T}
    color::Symbol
    name::String
end
function VirtualGate(sp::AbstractVector, ep::AbstractVector,
    n::AbstractVector, color::Symbol, name::String = "Virtual Gate")
    T = eltype(sp) <: Integer ? Float64 : eltype(sp)
    n = normalize(n)
    d = dot(n, ep-sp)
    if abs(d) > 10eps(T)
        error("Normal vector is not actually normal to the wall: dot = $d")
    end
    w = norm(ep - sp)
    center = @. (ep+sp)/2
    return VirtualGate{T}(SVector{2,T}(sp), SVector{2,T}(ep), SVector{2,T}(n), w, SVector{2,T}(center), color, name)
end

import DynamicalBilliards.specular!
@inline function specular!(p::AbstractParticle{T}, r::VirtualGate{T})::Void where {T}
    n = normalvec(r, p.pos)
    φ = atan2(n[2], n[1]) + 0.95(π*rand() - π/2) #this cannot be exactly π/2
    p.vel = SVector{2,T}(cos(φ), sin(φ))
    return nothing
end

import DynamicalBilliards.translation
translation(w::VirtualGate, vec) = VirtualGate(w.sp + vec, w.ep + vec, w.normal, w.color)

import DynamicalBilliards.collisiontime
function collisiontime(p::Particle{T}, w::VirtualGate{T})::T where {T}
    n = normalvec(w, p.pos)
    denom = dot(p.vel, n)
    # case of velocity pointing away of wall:
    denom ≥ 0 && return Inf
    posdot = dot(w.sp-p.pos, n)
    # Case of particle starting behind finite wall:
    posdot ≥ 0 && return Inf
    colt = posdot/denom
    intersection = p.pos .+ colt .* p.vel
    dfc = norm(intersection - w.center)
    if dfc > w.width/2
        return Inf
    else
        return colt
    end
end

@inline function project_to_line(point, c, n)
    posdot = dot(c - point, n)
    intersection = point + posdot*n
end

import DynamicalBilliards.distance_init
function distance_init(pos::SVector{2,T}, w::VirtualGate{T})::T where {T}
    n = normalvec(w, pos)
    posdot = dot(w.sp - pos, n)
    if posdot ≥ 0 # I am behind wall
        intersection = project_to_line(pos, w.center, n)
        dfc = norm(intersection - w.center)
        if dfc > w.width/2
            return +1.0 # but not directly behind
        else
            return -1.0
        end
    end
    v1 = pos - w.sp
    dot(v1, n)
end

max_stars_repo_name: gbarsih/Representations-in-Robotics

max_issues_repo_name: gbarsih/Representations-in-Robotics

max_forks_repo_name: gbarsih/Representations-in-Robotics


text:
#
# Load DFTK-compatible structural information from an external file
# Relies on ASE and other external libraries to do the parsing
#
using PyCall

function load_from_file(payload::Function, file::AbstractString)
    ase = pyimport_e("ase")
    ispynull(ase) && error("Install ASE to load data from exteral files")
    payload(pyimport("ase.io").read(file))
end
load_lattice(file::AbstractString)   = load_from_file(load_lattice,   file)
load_atoms(file::AbstractString)     = load_from_file(load_atoms,     file)
load_positions(file::AbstractString) = load_from_file(load_positions, file)
load_magnetic_moments(file::AbstractString) = load_from_file(load_magnetic_moments, file)

function load_lattice(pyobj::PyObject)
    mg = pyimport_e("pymatgen")
    ase = pyimport_e("ase")

    if !ispynull(mg)
        Lattice   = pyimport("pymatgen.core.lattice").Lattice
        Structure = pyimport("pymatgen.core.structure").Structure
        if any(pyisinstance(pyobj, s) for s in (Lattice, Structure))
            return load_lattice_pymatgen(pyobj)
        end
    end

    if !ispynull(ase)
        ase_supported = (ase.Atoms, ase.cell.Cell)
        if any(pyisinstance(pyobj, s) for s in ase_supported)
            return load_lattice_ase(pyobj)
        end
    end

    error("load_lattice not implemented for python type $pyobj")
end

function load_atoms(pyobj::PyObject)
    mg = pyimport_e("pymatgen")
    ase = pyimport_e("ase")

    if !ispynull(mg)
        if pyisinstance(pyobj, pyimport("pymatgen.core.structure").Structure)
            return load_atoms_pymatgen(pyobj)
        end
    end
    if !ispynull(ase) && pyisinstance(pyobj, ase.Atoms)
        return load_atoms_ase(pyobj)
    end

    error("load_atoms not implemented for python type $pyobj")
end

function load_positions(pyobj::PyObject)
    mg = pyimport_e("pymatgen")
    ase = pyimport_e("ase")

    if !ispynull(mg)
        if pyisinstance(pyobj, pyimport("pymatgen.core.structure").Structure)
            return load_positions_pymatgen(pyobj)
        end
    end
    if !ispynull(ase) && pyisinstance(pyobj, ase.Atoms)
        return load_positions_ase(pyobj)
    end

    error("load_positions not implemented for python type $pyobj")
end

function load_magnetic_moments(pyobj::PyObject)
    ase = pyimport_e("ase")
    if !ispynull(ase) && pyisinstance(pyobj, ase.Atoms)
        return load_magnetic_moments_ase(pyobj)
    end
    error("load_magnetic_moments not implemented for python type $pyobj")
end


text:
export export_network

getLabel(node::SumNode) = "\"+\""
getLabel(node::ProductNode) = "<×>"
getLabel(node::Leaf) = "\"X$(scope(node))\""
getLabel(node::IndicatorNode) = "\"X$(scope(node)) = $(node.value)\""

getShape(node::SPNNode) = "circle"
getShape(node::IndicatorNode) = "doublecircle"

getFontSize(node::SPNNode) = 40
getFontSize(node::IndicatorNode) = 14

getParameters(node::SPNNode) = "nothing"
getParameters(node::IndicatorNode) = "\"value=$(node.value)\""

function isDegenerated(node::Node, nodeObsDegeneration)
    return nodeObsDegeneration ? !(hasscope(node) && hasobs(node)) : !hasscope(node)
end

function isDegenerated(node::Leaf, nodeObsDegeneration)
    return all(map(p -> isDegenerated(p, nodeObsDegeneration), parents(node)))
end

function export_network(spn::SumProductNetwork,
                       filename::AbstractString;
                       nodeObsDegeneration = false,
                       excludeDegenerated = false)

    gstring = ["digraph SPN {"]
    push!(gstring, "node [margin=0, width=0.7, fixedsize=true];")

    if excludeDegenerated
        @warn("Excluding degenerated nodes, these nodes cannot be recovered from the DOT file!")
    end

    for node in values(spn)
        label = getLabel(node)
        shape = getShape(node)
        fontsize = getFontSize(node)
        parameters = getParameters(node)
        nodeType = "\"$(string(typeof(node)))\""
        nodeScope = "\"$(string(scope(node)))\""
        nodeNumObs = isa(node, Node) ? string(sum(node.obsVec)) : "\"\""

        degeneratedNode = isDegenerated(node, nodeObsDegeneration)
        style = !degeneratedNode ? "" : "style=dotted,"

        if excludeDegenerated && degeneratedNode
            continue
        end

        nodestring = "$(spn.idx[node.id]) [label=$(label), shape=$(shape), $(style) fontsize=$(fontsize), "*
                        "nodeType=$(nodeType), nodeScope=$(nodeScope), nodeNumObs=$(nodeNumObs), "*
                        "nodeParameters=$(parameters)];"

        push!(gstring, nodestring)
    end

    for node in filter(n -> isa(n, Node), values(spn))
        for (ci, child) in enumerate(children(node))
            degeneratedChild = isDegenerated(child, nodeObsDegeneration)
            degeneratedNode = isDegenerated(node, nodeObsDegeneration)
            style = !degeneratedChild ? "" : "style=dotted,"

            if excludeDegenerated & (degeneratedNode || degeneratedChild)
                continue
            end

            if hasweights(node)
                w = weights(node)[ci]
                #logw = @sprintf("%0.4f", logw)
                w = @sprintf("%0.2f", w)
                push!(gstring, "$(spn.idx[node.id]) -> $(spn.idx[child.id]) [label=$(w), $(style) logWeight=\"$(logweights(node))\"];")

            else
                style = !degeneratedChild ? "" : "style=dotted"
                push!(gstring, "$(spn.idx[node.id]) -> $(spn.idx[child.id]) [$(style)];")
            end
        end
    end

    push!(gstring, "}")

    open(filename, "w") do f
        write(f, join(gstring, '\n'))
    end
end


text:
"""
    Constellation

Constellation type for simulating digital waveforms.
"""
struct Constellation
    label   :: String
    symbols :: Vector{ComplexF64}  # symbol alphabet
end

avgenergy(v) = mean(abs2.(v))
avgenergy(c::Constellation) = avgenergy(c.symbols)
Base.length(c::Constellation) = length(c.symbols)
Base.getindex(c::Constellation, index) = c.symbols[index]
Base.rand(c::Constellation, n::Integer) = rand(c.symbols, n)
Base.rand(c::Constellation) = rand(c.symbols, 1)
#Base.show(c::Constellation) = plot(c.symbols, w="p")

"""
    pam(; M = 2, Es = 1.0, ϕ = 0.0) -> Constellation

Return a PAM constellation with `M` elements and average symbol energy `Es`.
The constellation is rotated ϕ radians.
"""
function pam(; M = 2, Es = 1.0, ϕ = 0.0)
    isinteger(log2(M)) || error("M must be a power of 2 (M = $M)")
    v = M/2
    base = range(-(2v-1), (2v-1), step = 2)
    E = avgenergy(base)
    symbols = exp(1im*ϕ)*sqrt(Es/E).*base
    return Constellation("$M-PAM", collect(symbols))
end

function squareqam(; M = 4, Es = 1.0, ϕ = 0.0)
    iseven(log2(M)) || error("M must be an even power of 2 (M = $M)")
    M2 = Int(log2(M))
    c = real(pam(M=M2).symbols)
    symbols = [ComplexF64(i, q) for i in c for q in c]
    E = avgenergy(symbols)
    symbols .*= sqrt(Es/E)*exp(1im*ϕ)
    return Constellation("$M-QAM", symbols)
end

"""
    qam(; M = 4, Es = 1.0, ϕ = 0.0)

Return a QAM constellation with `M` elements and average symbol energy `Es`.
The constellation is rotated ϕ radians.

Note: only square constellations (where `M = 2^(2k)``) are supported.
"""
function qam(; M = 4, kwargs...)
    if iseven(log2(M))
        return squareqam(; M = M, kwargs...)
    else
        error("Only square QAM constellations are supported (M = $M)")
    end
end

"""
    psk(; M = 4, Es = 1.0, ϕ = 0.0)

Return a PSK constellation with `M` elements and average symbol energy `Es`.
The constellation is rotated ϕ radians.
"""
function psk(; M = 4, Es = 1.0, ϕ = 0.0)
    isinteger(log2(M)) || error("M must be a power of 2 (M = $M)")
    θ = range(0, step = 2π/M, length = M)
    symbols = exp.(1im*θ)
    E = avgenergy(symbols)
    symbols .*= sqrt(Es/E)*exp(1im*ϕ)
    return Constellation("$M-PSK", symbols)
end


text:
var_desc = Dict(
    "TEMP" => Dict(
        "long_name" => "Temperature",
        "units"     => "degC",
    ),

    "SALT" => Dict(
        "long_name" => "Salinity",
        "units"     => "kg / m^3 (PSU)",
    ),

)

max_stars_repo_name: meteorologytoday/IOM

max_issues_repo_name: meteorologytoday/IOM

max_forks_repo_name: meteorologytoday/IOM


text:
struct LorentzVectorCyl{T <: AbstractFloat}
    pt::T
    eta::T
    phi::T
    mass::T
end

Base.show(io::IO, v::LorentzVectorCyl) = print(io, "$(typeof(v))(pt=$(v.pt), eta=$(v.eta), phi=$(v.phi), mass=$(v.mass))")

Base.broadcastable(v::LorentzVectorCyl) = Ref(v)

"""
    zero(LorentzVectorCyl)
    zero(LorentzVectorCyl{T})
Constructs a zero four-vector.
"""
Base.zero(::Type{LorentzVectorCyl{T}}) where T = LorentzVectorCyl{T}(zero(T), zero(T), zero(T), zero(T))
Base.zero(::Type{LorentzVectorCyl}) = zero(LorentzVectorCyl{Float64})

px(v::LorentzVectorCyl) = v.pt * cos(v.phi)
py(v::LorentzVectorCyl) = v.pt * sin(v.phi)
pz(v::LorentzVectorCyl) = v.pt * sinh(v.eta)
energy(v::LorentzVectorCyl) = sqrt(px(v)^2 + py(v)^2 + pz(v)^2 + v.mass^2)

function *(v::LorentzVectorCyl{T}, k::Real) where T
    LorentzVectorCyl{T}(v.pt*k, v.eta, v.phi, v.mass*k)
end

function +(v1::LorentzVectorCyl{T}, v2::LorentzVectorCyl{W}) where {T,W}
    m1, m2 = max(v1.mass, zero(v1.pt)), max(v2.mass, zero(v2.pt))

    px1, px2 = px(v1), px(v2)
    py1, py2 = py(v1), py(v2)
    pz1, pz2 = pz(v1), pz(v2)
    e1 = sqrt(px1^2 + py1^2 + pz1^2 + m1^2)
    e2 = sqrt(px2^2 + py2^2 + pz2^2 + m2^2)

    sumpx = px1+px2
    sumpy = py1+py2
    sumpz = pz1+pz2

    ptsq = sumpx^2 + sumpy^2
    pt = sqrt(ptsq)
    eta = asinh(sumpz/pt)
    phi = atan(sumpy, sumpx)
    mass = sqrt(max(muladd(m1, m1, m2^2) + 2*e1*e2 - 2*(muladd(px1, px2, py1*py2) + pz1*pz2), zero(v1.pt)))
    return LorentzVectorCyl(pt,eta,phi,mass)
end

function fast_mass(v1::LorentzVectorCyl, v2::LorentzVectorCyl)
    # Calculate mass directly. Same as (v1+v2).mass except
    # this skips the intermediate pt, eta, phi calculations.
    # ~4x faster than (v1+v2).mass
    pt1, pt2 = v1.pt, v2.pt
    eta1, eta2 = v1.eta, v2.eta
    phi1, phi2 = v1.phi, v2.phi
    m1, m2 = v1.mass, v2.mass

    # note, massless approximation is
    # mass = sqrt(max(2*pt1*pt2*(cosh(eta1-eta2) - cos(phi1-phi2)), zero(pt1)))

    # @fastmath sinh() is very slow in some cases
    sinheta1 = sinh(eta1)
    sinheta2 = sinh(eta2)
    tpt12 = 2*pt1*pt2
    return @fastmath sqrt(max(muladd(m1, m1, m2^2)
        + 2*sqrt((pt1^2*(1+sinheta1^2) + m1^2)*(pt2^2*(1+sinheta2^2) + m2^2))
        - tpt12*sinheta1*sinheta2
        - tpt12*cos(phi1-phi2), zero(pt1)))
end

# https://root.cern.ch/doc/v606/GenVector_2VectorUtil_8h_source.html#l00061
@inline function deltaphi(v1::LorentzVectorCyl, v2::LorentzVectorCyl)
    dphi = v2.phi - v1.phi
    dphi = dphi - 2*pi*(dphi > pi) + 2*pi*(dphi <= -pi)
    return dphi
end
const Δϕ = deltaphi

@inline function deltaeta(v1::LorentzVectorCyl, v2::LorentzVectorCyl)
    return v2.eta - v1.eta
end
const Δη = deltaeta

@inline function deltar2(v1::LorentzVectorCyl, v2::LorentzVectorCyl)
    dphi = deltaphi(v1,v2)
    deta = deltaeta(v1,v2)
    return muladd(dphi, dphi, deta^2)
end
@inline deltar(v1::LorentzVectorCyl, v2::LorentzVectorCyl) = sqrt(deltar2(v1, v2))
const ΔR = deltar

function tocartesian(v::LorentzVectorCyl)
    x = px(v)
    y = py(v)
    z = pz(v)
    t = energy(v)
    (t=t, x=x, y=y, z=z)
end

function fromcartesian(t, x, y, z)
    pt = sqrt(muladd(x, x, y^2))
    eta = asinh(z/pt)
    phi = atan(y, x)
    mass = sqrt(max(t^2 - x^2 - y^2 - z^2, 0))
    return LorentzVectorCyl(pt,eta,phi,mass)
end


text:
module AbstractBijections
  using Dictionaries

  import Base: insert!, getindex, inv, length, show, copy

  export Bijection, AbstractBijection, bijection, domain, image, domain_eltype, image_eltype

  abstract type AbstractBijection{D,I} end

  _not_implemented() = error("Not implemented")

  # Apply the Bijective function to an element
  # `x ∈ domain(f)`, resulting in an element `y ∈ image(f)`.
  apply(f::AbstractBijection, x) = _not_implemented()
  inv(f::AbstractBijection) = _not_implemented()
  domain(f::AbstractBijection) = _not_implemented()
  image(f::AbstractBijection) = _not_implemented()
  insert!(f::AbstractBijection, x, y) = _not_implemented()

  apply_inv(f::AbstractBijection, y) = apply(inv(f), y)

  getindex(f::AbstractBijection, x) = apply(f, x)
  length(f::AbstractBijection) = length(domain(f))
  domain_eltype(f::AbstractBijection) = eltype(domain(f))

  function show(io::IO, mime::MIME"text/plain", f::AbstractBijection)
    println(io, typeof(f))
    for (x, y) in zip(domain(f), image(f))
      show(io, x)
      print(io, " ↔ ")
      show(io, y)
      println(io)
    end
    return nothing
  end

  show(io::IO, f::AbstractBijection) = show(io, MIME"text/plain"(), f)

  struct Bijection{D,I} <: AbstractBijection{D,I}
    f::Dictionary{D,I}
    finv::Dictionary{I,D}
  end

  copy(f::Bijection) = Bijection(copy(f.f), copy(f.finv))

  function Bijection(domain, image)
    f = Dictionary(domain, image)
    finv = Dictionary(image, domain)
    return Bijection(f, finv)
  end
  Bijection(image) = Bijection(eachindex(image), image)

  function bijection(domain_image)
    f = dictionary(domain_image)
    finv = Dictionary(f.values, f.indices)
    return Bijection(f, finv)
  end

  apply(f::Bijection, x) = f.f[x]
  inv(f::Bijection) = Bijection(f.finv, f.f)

  domain(f::Bijection) = f.finv.values
  image(f::Bijection) = f.f.values

  function insert!(f::Bijection, x, y)
    insert!(f.f, x, y)
    insert!(f.finv, y, x)
    return f
  end
end # module AbstractBijections


text:
#######
# stock names is a vector of all the stock names corresponding to columns in R
# p is a matrix of all the prices
# R is the return matrix
# dates is an array of all the dates corresponding to specific rows in R

### Loads in the matrix
using DataArrays, DataFrames

prices_data = readtable("prices.dat")
Symbols = Base.names(prices_data);

#Extract the stock symbols from the data frame
stockNames = String[];
for symb=Symbols
	strSymb = string(symb)
	if strSymb != "Row" && strSymb != "Dates"
		push!(stockNames, strSymb)
	end
end
stock_symbols = map(Symbol, stockNames)

#Extract the R matrix and dates
p = convert(Array,prices_data[stock_symbols]);
dates = prices_data[Symbol("Dates")] ;
R = diff(p) ./ p[1:end-1, :];
T,n = size(R);


text:
module Factorial

f(x::Int64) = x <= 1 ? 1 : x * f(x - 1)

end

@ctx (true, true, true) struct RecursionMix end
allow(ctx::RecursionMix, m::Module) = m == Factorial
debug(::RecursionMix) = false

swap(e) = e
function swap(e::Expr)
    new = MacroTools.postwalk(e) do s
        s isa Expr || return s
        s.head == :call || return s
        s.args[1] == Base.:(*) || return s
        return Expr(:call, Base.:(+), e.args[2:end]...)
    end
    return new
end

function transform(::RecursionMix, b)
    for (v, st) in b
        replace!(b, v, swap(st))
    end
    return b
end

@testset "Recursion" begin
    Mixtape.@load_call_interface()
    @test call(RecursionMix(), Factorial.f, 10) == 55
end

using Core.Compiler: NativeInterpreter
display(Mixtape.CACHE[NativeInterpreter])


text:
"""
find_gaps(time, Δt = median(diff(time)); radius = 0.1Δt)

Find gaps in a time vector and encode them into a matrix S. S' * y results in the original signal with zero filling.

# Arguments

* time:     A vector that encodes the time (eg. seconds, samples, ...)
* Δt:       The expected distance between two measurements, defaults to median(diff(time))
* radius:   If the signal is not non-uniformly sampled, the radius can be used to search in a broader area around the expected sample.
"""
function find_gaps(time, Δt = median(diff(time)); radius = 0.1Δt)
    present = Bool[]
    diffTime = diff(time)
    currentGapLen = 1
    push!(present, true)
    for difft in diffTime
        while currentGapLen * Δt + radius < difft
            push!(present, false)
            currentGapLen += 1
        end
        push!(present, true)
        currentGapLen = 1
    end
    S = I(length(present))[present, :]
    return S
end

"""
simulate_data(f; percentage_gaps = 0.2, max_gap_ratio = 0.3, percentage_outliers = 0.05, σ_noise = 0.2, outlier_strength = 1)

Put some gaps and outlier in the smooth signal f. The underlying model takes the form:
Sy = Sf + x + w

# Keyword Arguments

* percentage_gaps: The percentage of gaps ∈ [0,1]
* max_gap_ratio: Maximal allowed ratio between gap length to total gap length ∈ ]0,1]
* percentage_outliers: The percentage of outliers ∈ [0,1]
* outlier_strength: Strength of the gaussian outlier distribution
* σ_noise: σ of the additive gaussian noise w

# Returns

(y, x, w, S)
"""
function simulate_data(f; percentage_gaps = 0.4, max_gap_ratio = 0.1, percentage_outliers = 0.05, σ_noise = 0.2, outlier_strength = 1)
    N = length(f)
    # Gaps
    present = place_gaps(N, percentage_gaps = percentage_gaps, max_gap_ratio = max_gap_ratio)
    S = I(N)[present, :]
    # Outliers
    x = zeros(N)
    outliers = sample(present, round(Int, percentage_outliers * N), replace = false)
    x[outliers] += rand(Normal(0, outlier_strength), length(outliers))
    # Noise
    w = rand(Normal(0, σ_noise), N)
    y = f .+ x .+ w
    return (y, x, w, S)
end

"""
Get the indices where data is present from the encoding matrix S
"""
get_present(S) = [1:size(S,2);][diag(S' * S) .== 1]

# Distribute the gaps onto the signal, return which samples are present
function place_gaps(N; percentage_gaps = 0.2, max_gap_ratio = 0.3)
    present = collect(1:N)
    total_gaps = round(Int, percentage_gaps * N)
    max_gap = max_gap_ratio * total_gaps
    while total_gaps > 0
        gap_len = sample(1:min(max_gap, total_gaps))
        start = sample(1:(length(present) - gap_len))
        filter!(x -> !(x in start:(start + gap_len - 1)), present)
        total_gaps = length(present) - (N - round(Int, percentage_gaps * N))
    end
    return present
end

# L1 norm
norm1(x) = sum(abs, x)

# L2 norm
norm2(x) = sqrt(sum(abs2, x))


text:
include("../anderson.jl")
savefig(plot_anderson(@__DIR__), "damping.svg")