text:
# Copyright (c) 2018: Matthew Wilhelm & Matthew Stuber.
# This code is licensed under MIT license (see LICENSE.md for full details)
#############################################################################
# EAGO
# A development environment for robust and global optimization
# See https://github.com/PSORLab/EAGO.jl
#############################################################################
# src/eago_semiinfinite/semi_infinite.jl
# Defines high level interface for SIP routines.
#############################################################################

include("types.jl")
include("subproblems.jl")

function sip_solve(t::ExtensionType, alg::AbstractSIPAlgo, subresult, prob, result, cb)
    error("Algorithm $t not supported by explicit_sip_solve.")
end

"""
    sip_solve

Solve an SIP with decision variable bounds `x_l` to `x_u`, uncertain variable
bounds `p_l` to `p_u`, an objective function of `f`, and `gSIP` seminfiniite
constraint(s).
"""
function sip_solve(alg::T, x_l::Vector{Float64}, x_u::Vector{Float64},
                   p_l::Vector{Float64}, p_u::Vector{Float64},
                   f::Function, gSIP::Vector{Any}; d = DefaultExt(), kwargs...) where {T <: AbstractSIPAlgo}

    @assert length(p_l) == length(p_u)
    @assert length(x_l) == length(x_u)

    prob = SIPProblem(x_l, x_u, p_l, p_u, gSIP, kwargs)
    subresult = SIPSubResult(prob.nx, prob.np, prob.nSIP, prob.abs_tolerance)
    result = SIPResult(prob.nx, prob.np)
    cb = SIPCallback(f, gSIP)

    sip_solve!(d, alg, subresult, prob, result, cb)
    return result
end

function sip_solve(alg::T, x_l::Vector{Float64}, x_u::Vector{Float64},
                   p_l::Vector{Float64}, p_u::Vector{Float64},
                   f::Function, gSIP::Function; d = DefaultExt(), kwargs...) where {T <: AbstractSIPAlgo}
    return sip_solve(alg, x_l, x_u, p_l, p_u, f, [gSIP], d = d, kwargs...)
end

include("nonconvex_algorithms/sip_hybrid.jl")
include("nonconvex_algorithms/sip_res_rev.jl")
include("nonconvex_algorithms/sip_res.jl")

ext: jl


text:
# Given a DNA strand, its transcribed RNA strand is formed by replacing each nucleotide with its complement:
# G -> C, C -> G, T -> A, A -> U
function to_rna(dna::AbstractString)
    typeof(match(r"^[GCTA]+$", dna)) == Nothing && error("Invalid RNA strand")
    # Define character associations
    a_nucleotides = Dict('G'=>'C', 'C'=>'G', 'T'=>'A', 'A'=>'U')
    # Replace characters using dictionary
    map((x)->a_nucleotides[x], dna)
end


text:
function plot_projection2D!(p, proj; l = 10, spkwargs...)

    def_spkwargs = (;label = "", alpha = 0.5)

    flx1s = sort(collect(keys(proj)))
    rflx1s = reverse(flx1s)
    flx2lbs = [first(proj[flx1]) for flx1 in flx1s]
    flx2ubs = [last(proj[flx1]) for flx1 in rflx1s]

    s = Shape([flx1s; rflx1s], [flx2lbs; flx2ubs])
    plot!(p, [s]; label = "", def_spkwargs..., spkwargs...)

    return p
end

## ----------------------------------------------------------------------------
function plot_projection2D!(p::AbstractPlot, model::ChU.MetNet,
        ider1::ChU.IDER_TYPE, ider2::ChU.IDER_TYPE;
        l = 10, spkwargs...
    )

    proj = ChLP.projection2D(model, ider1, ider2; l)
    plot_projection2D!(p, proj, l, spkwargs...)

end


text:
using BinaryBuilder

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

sources_linux_x64 = [
    ArchiveSource("https://developer.nvidia.com/compute/machine-learning/cudnn/secure/8.0.3.33/10.2_20200825/cudnn-10.2-linux-x64-v8.0.3.33.tgz",
                  "b3d487c621e24b5711983b89bb8ad34f0378bdbf8a1a4b86eefaa23b19956dcc")
]
sources_linux_ppc64le = [
    ArchiveSource("https://developer.nvidia.com/compute/machine-learning/cudnn/secure/8.0.3.33/10.2_20200825/cudnn-10.2-linux-ppc64le-v8.0.3.33.tgz",
                  "ff22c9c37af191c9104989d784427cde744cdde879bfebf3e4e55ca6a9634a11")
]
sources_windows = [
    ArchiveSource("https://developer.nvidia.com/compute/machine-learning/cudnn/secure/8.0.3.33/10.2_20200825/cudnn-10.2-windows10-x64-v8.0.3.33.zip",
                  "45ffcb6ef3a493995d442a8503186f54c5147a1a3c94fdf9e40846a2ce53b64b")
]

include("../common.jl")

max_stars_repo_name: mortenpi/Yggdrasil

max_issues_repo_name: mortenpi/Yggdrasil

max_forks_repo_name: mortenpi/Yggdrasil


text:
function movetoside!(r,side)
    count=0
    while isborder(r,side)==0
        move!(r,side)
        count=count+1
    end
    return count
end

function corner!(r)
home=[]
    while isborder(r,Nord)==0 || isborder(r,West)==0
        x=movetoside!(r,Nord)
        y=movetoside!(r,West)
        push!(home,x)
        push!(home,y)
    end
    return home
end

function countmove!(r,side,count)
    for i in 1:count
        if isborder(r,side)==0
        move!(r,side)
        end
    end
end

function gohome!(r,home)
    reverse!(home)
    for i = 1:length(home)
        if i%2==0
            countmove!(r,Sud,home[i])
        else
            countmove!(r,Ost,home[i])
        end
    end
end
function Points!(r)
    while isborder(r,Sud)==0
        while isborder(r,Ost)==0
         putmarker!(r)
          move!(r,Ost)
        end

    if isborder(r,Ost)==1
        putmarker!(r)
        if isborder(r,Sud)==0
        move!(r,Sud)
        else
        while isborder(r,West)==0
            move!(r,West)
        end
        end
        while isborder(r,West)==0
        putmarker!(r)
        move!(r,West)
        end
    end
    putmarker!(r)
    if isborder(r,Sud)==0
    move!(r,Sud)

    end
if isborder(r,Sud)==1
    while isborder(r,Ost)==0
        putmarker!(r)
         move!(r,Ost)
       end
       putmarker!(r)
end
end
putmarker!(r)

while isborder(r,West)==0
    move!(r,West)
end
while isborder(r,Nord)==0
    move!(r,Nord)
end
end

function task_30(r)
home=[]
home=corner!(r)
Points!(r)
gohome!(r,home)
end

task_30(r)


text:
"""
Multithreaded CSR Matrix-Vector Product
"""
module MtSpMV

using SparseArrays, LinearAlgebra

using Base: promote_op

import Base: *
import LinearAlgebra: mul!, matprod

export mul!, *

include("parallellinalg.jl")

end


text:
"""
    Libcint

Minimal wrap around the integral library libcint. This module exposes
libcint functions to the Julia interface.
"""
module Libcint

export cint1e_kin_sph!, cint1e_nuc_sph!, cint1e_ovlp_sph!, cint2c2e_sph!, cint2e_sph!, cint3c2e_sph!
export cint1e_ipkin_sph!, cint1e_ipnuc_sph!, cint1e_ipovlp_sph!, cint2e_ip1_sph!, cint1e_r_sph!

const LIBCINT = joinpath(@__DIR__, "../deps/libcint")

function CINTcgtos_spheric(bas_id, bas)
    @ccall LIBCINT.CINTcgtos_spheric(bas_id::Cint, bas::Ptr{Cint})::Cint
end

function cint1e_ovlp_sph!(buf, shls, atm, natm, bas, nbas, env)
    @ccall LIBCINT.cint1e_ovlp_sph(
                                    buf  :: Ptr{Cdouble},
                                    shls :: Ptr{Cint},
                                    atm  :: Ptr{Cint},
                                    natm :: Cint,
                                    bas  :: Ptr{Cint},
                                    nbas :: Cint,
                                    env  :: Ptr{Cdouble}
                                   )::Cvoid
end

function cint1e_kin_sph!(buf, shls, atm, natm, bas, nbas, env)
    @ccall LIBCINT.cint1e_kin_sph(
                                    buf  :: Ptr{Cdouble},
                                    shls :: Ptr{Cint},
                                    atm  :: Ptr{Cint},
                                    natm :: Cint,
                                    bas  :: Ptr{Cint},
                                    nbas :: Cint,
                                    env  :: Ptr{Cdouble}
                                   )::Cvoid
end

function cint1e_nuc_sph!(buf, shls, atm, natm, bas, nbas, env)
    @ccall LIBCINT.cint1e_nuc_sph(
                                    buf  :: Ptr{Cdouble},
                                    shls :: Ptr{Cint},
                                    atm  :: Ptr{Cint},
                                    natm :: Cint,
                                    bas  :: Ptr{Cint},
                                    nbas :: Cint,
                                    env  :: Ptr{Cdouble}
                                   )::Cvoid
end

function cint2e_sph!(buf, shls, atm, natm, bas, nbas, env)
    opt = Ptr{UInt8}(C_NULL)
    @ccall LIBCINT.cint2e_sph(
                                    buf  :: Ptr{Cdouble},
                                    shls :: Ptr{Cint},
                                    atm  :: Ptr{Cint},
                                    natm :: Cint,
                                    bas  :: Ptr{Cint},
                                    nbas :: Cint,
                                    env  :: Ptr{Cdouble},
                                    opt :: Ptr{UInt8},
                                   )::Cvoid
end

function cint2c2e_sph!(buf, shls, atm, natm, bas, nbas, env)
    opt = Ptr{UInt8}(C_NULL)
    @ccall LIBCINT.cint2c2e_sph(
                                    buf  :: Ptr{Cdouble},
                                    shls :: Ptr{Cint},
                                    atm  :: Ptr{Cint},
                                    natm :: Cint,
                                    bas  :: Ptr{Cint},
                                    nbas :: Cint,
                                    env  :: Ptr{Cdouble},
                                    opt :: Ptr{UInt8},
                                   )::Cvoid
end

function cint3c2e_sph!(buf, shls, atm, natm, bas, nbas, env)
    opt = Ptr{UInt8}(C_NULL)
    @ccall LIBCINT.cint3c2e_sph(
                                    buf  :: Ptr{Cdouble},
                                    shls :: Ptr{Cint},
                                    atm  :: Ptr{Cint},
                                    natm :: Cint,
                                    bas  :: Ptr{Cint},
                                    nbas :: Cint,
                                    env  :: Ptr{Cdouble},
                                    opt :: Ptr{UInt8},
                                   )::Cvoid
end

function cint1e_ipovlp_sph!(buf, shls, atm, natm, bas, nbas, env)
    @ccall LIBCINT.cint1e_ipovlp_sph(
                                    buf  :: Ptr{Cdouble},
                                    shls :: Ptr{Cint},
                                    atm  :: Ptr{Cint},
                                    natm :: Cint,
                                    bas  :: Ptr{Cint},
                                    nbas :: Cint,
                                    env  :: Ptr{Cdouble}
                                   )::Cvoid
end

function cint1e_ipkin_sph!(buf, shls, atm, natm, bas, nbas, env)
    @ccall LIBCINT.cint1e_ipkin_sph(
                                    buf  :: Ptr{Cdouble},
                                    shls :: Ptr{Cint},
                                    atm  :: Ptr{Cint},
                                    natm :: Cint,
                                    bas  :: Ptr{Cint},
                                    nbas :: Cint,
                                    env  :: Ptr{Cdouble}
                                   )::Cvoid
end

function cint1e_ipnuc_sph!(buf, shls, atm, natm, bas, nbas, env)
    @ccall LIBCINT.cint1e_ipnuc_sph(
                                    buf  :: Ptr{Cdouble},
                                    shls :: Ptr{Cint},
                                    atm  :: Ptr{Cint},
                                    natm :: Cint,
                                    bas  :: Ptr{Cint},
                                    nbas :: Cint,
                                    env  :: Ptr{Cdouble}
                                   )::Cvoid
end

function cint2e_ip1_sph!(buf, shls, atm, natm, bas, nbas, env)
    opt = Ptr{UInt8}(C_NULL)
    @ccall LIBCINT.cint2e_ip1_sph(
                                    buf  :: Ptr{Cdouble},
                                    shls :: Ptr{Cint},
                                    atm  :: Ptr{Cint},
                                    natm :: Cint,
                                    bas  :: Ptr{Cint},
                                    nbas :: Cint,
                                    env  :: Ptr{Cdouble},
                                    opt :: Ptr{UInt8},
                                   )::Cvoid
end

function cint1e_r_sph!(buf, shls, atm, natm, bas, nbas, env)
    @ccall LIBCINT.cint1e_r_sph(
                                    buf  :: Ptr{Cdouble},
                                    shls :: Ptr{Cint},
                                    atm  :: Ptr{Cint},
                                    natm :: Cint,
                                    bas  :: Ptr{Cint},
                                    nbas :: Cint,
                                    env  :: Ptr{Cdouble}
                                   )::Cvoid
end

end #module


text:
# Rijke tube model
# Reference: Huhn and Magri 2020
using LinearAlgebra
#using DiffEqSensitivity, OrdinaryDiffEq, Zygote
include("cheb.jl")
# s = [beta, tau]
c1 = 0.1
c2 = 0.06
xf = 0.2
Ng = 10
tNg = 2*Ng + 1
Nc = 10
N = 2*Ng + Nc
j = LinRange(1,Ng,Ng)
jpi = @. pi.*j
cjpixf = @. cos.(xf.*jpi)
sjpixf = @. sin.(xf.*jpi)
zetaj = @. c1.*j.*j .+ c2.*
		sqrt.(j)
coeffs = [-1.0, 0.0, 1.75e3, 6.2e-12, -7.5e6]
tau = 0.2
dt = 5.e-1/Nc*tau
function qfun(t)
	if abs(t + 1.0) > 0.005
		return sqrt(abs(1.0 + t)) - 1.0
	end
	q = 0.
	for i = 1:5
		q += coeffs[i]*((1+t)^(i-1))
	end
	return q
end

D = cheb_diff_matrix(Nc)
function Rijke(u0::Array{Float64,1},
			   s::Array{Float64,1}=[7.0, 0.2],
			   n::Int64=1)
	uDot1 = zeros(N)
	uDot2 = zeros(N)
	uDot3 = zeros(N)
	beta, tau = s
	#ufMean = zeros(n)
	u = copy(u0)
	for i = 1:n
		#ufMean[i] = dot(u0[1:Ng],cjpixf)
		un = copy(u)
		f!(uDot1, u, s, i)
		@. u = u + (1/2)*dt*uDot1
		f!(uDot2, u, s, i)
		@. u = u - dt*uDot1 +
			2*dt*uDot2
		f!(uDot3, u, s, i)
		@. u = un + dt*((1/6).*uDot1 +
					 (2/3)*uDot2 + (1/6)*uDot3)
	end
	return u
end
function Rijke_ODE(u0::Array{Float64,1},
			   s::Array{Float64,1}=[7.0, 0.2],
			   n::Int64=1)
	t = n*dt
	prob = ODEProblem(f!, u0, (0.,t), s)
	sol = Array(solve(prob, Tsit5(),saveat=dt))
	return sol[:,end]
end
function f!(uDot, u, s, t)
	eta = view(u, 1:Ng)
	mu = view(u, Ng+1:2*Ng)
	v = view(u, tNg:N)
	beta, tau = s
	heat = beta*qfun(v[1])
	@. uDot[1:Ng] = jpi*mu
	@. uDot[Ng+1:2*Ng] = -jpi*eta - zetaj*mu - 2.0*heat*
						sjpixf
	uDot[tNg:N] .= -2.0/tau.*(D*[v;
					dot(eta, cjpixf)])[1:end-1]
end
function dfdbeta(u,s)
	n, d = size(u)
	beta, tau = s
	v1 = u[:,2*Ng + 1]
	dheat = qfun.(v1)
	return reshape([zeros(Ng*n);
					kron(-2.0.*ones(n).*dheat, sjpixf);
					zeros(Nc*n)], d, n)
end
function perturbation(u, s, eps=1.e-6)
	d, = size(u)
	beta, tau = s
	return (Rijke(u, [beta + eps, tau], 1) -
			Rijke(u, [beta - eps, tau], 1))/(2*eps)
end
function tau_perturbation(u, s, eps=1.e-6)
	d, = size(u)
	beta, tau = s
	return (Rijke(u, [beta, tau + eps], 1) -
			Rijke(u, [beta, tau - eps], 1))/(2*eps)
end

function dRijke(u, s, eps)
	d,  = size(u)
	dTu = zeros(d,d)
	u_p = zeros(d)
	u_m = zeros(d)
	v0 = zeros(d)
	for j = 1:d
		v0 .= zeros(d)
		v0[j] = 1.0
		u_p .= u .+ eps*v0
		u_m .= u .- eps*v0

		dTu[:,j] = (Rijke_ODE(u_p, s, 1) -
					Rijke_ODE(u_m, s, 1))/(2*eps)
	end
	return dTu
end
function dRijke_AD(u, s, eps=0.)
    du_ad = zeros(d,d)
    for j = 1:d
        v = zeros(d)
        du_ad[:,j] = Zygote.gradient(v->
                Rijke_ODE(u .+ v, s, 1)[j],v)[1]
    end
    return du_ad'
end

max_stars_repo_name: nishaChandramoorthy/energies

max_issues_repo_name: nishaChandramoorthy/energies

max_forks_repo_name: nishaChandramoorthy/energies


text:
import BenchPerf
import Random
using BenchmarkTools

function increment_at!(xs, indices)
    for i in indices
        @inbounds xs[i] += 1
    end
end

CACHE = Dict()
group = BenchmarkGroup()
for e in 5:22
    n = 2^e
    CACHE[e] = (xs = zeros(Int, n), indices = rand(1:n, n))
    group["n=$n"] = bench = @benchmarkable(
        increment_at!(xs, indices),
        setup = begin
            inputs = CACHE[$e]::$(typeof(CACHE[e]))
            xs = inputs.xs
            indices = inputs.indices
            fill!(xs, 0)
        end,
        # Larger `samples` (than the default: 10_000) to bound the benchmark by
        # the time limit.
        samples = 1_000_000,
    )
    # Manually tune the benchmark since `BenchmarkConfigSweeps` does not do it
    # ATM.
    if e < 15
        tune!(bench)
    end
end

include("../quick.jl")
maybe_quick!(group)

SUITE = BenchPerf.wrap(group; detailed = 1)


text:
using Gadfly
import Dolo
using AxisArrays

path = Dolo.pkg_path

model = Dolo.yaml_import(
        joinpath(path,
         "examples","models",
         "rbc_catastrophe.yaml")
)

Dolo.discretize(Dolo.MarkovChain, model.exogenous)

dp = Dolo.discretize(model.exogenous)

sol = Dolo.time_iteration(model, dp)

sim = Dolo.simulate(model, sol.dr, dp, T=100)
sim
plot(y=sim[Axis{:N}(1), Axis{:V}(:k)],Geom.line)

plot(y=sim[Axis{:N}(1), Axis{:V}(:mc_process)],Geom.line)

dp

Gadfly.plot(y=)


text:
#
# We use these definitions because when we have other SIMD operations with masks
# LLVM optimizes the masks better.
function truncate_mask!(instrs, input, W, suffix)
    mtyp_input = "i$(max(8,W))"
    mtyp_trunc = "i$(W)"
    str = if mtyp_input == mtyp_trunc
        "%mask.$(suffix) = bitcast $mtyp_input %$input to <$W x i1>"
    else
        "%masktrunc.$(suffix) = trunc $mtyp_input %$input to $mtyp_trunc\n%mask.$(suffix) = bitcast $mtyp_trunc %masktrunc.$(suffix) to <$W x i1>"
    end
    push!(instrs, str)
end
function zext_mask!(instrs, input, W, suffix)
    mtyp_input = "i$(max(8,W))"
    mtyp_trunc = "i$(W)"
    str = if mtyp_input == mtyp_trunc
        "%res.$(suffix) = bitcast <$W x i1> %$input to $mtyp_input"
    else
        "%restrunc.$(suffix) = bitcast <$W x i1> %$input to $mtyp_trunc\n%res.$(suffix) = zext $mtyp_trunc %restrunc.$(suffix) to $mtyp_input"
    end
    push!(instrs, str)
end
function binary_mask_op_instrs(W, op)
    mtyp_input = "i$(max(8,W))"
    instrs = String[]
    truncate_mask!(instrs, '0', W, 0)
    truncate_mask!(instrs, '1', W, 1)
    push!(instrs, "%combinedmask = $op <$W x i1> %mask.0, %mask.1")
    zext_mask!(instrs, "combinedmask", W, 1)
    push!(instrs, "ret $mtyp_input %res.1")
    join(instrs, "\n")
end
function binary_mask_op(W, U, op)
    instrs = binary_mask_op_instrs(W, op)
    quote
        $(Expr(:meta,:inline))
        Mask{$W}(llvmcall($instrs, $U, Tuple{$U, $U}, m1.u, m2.u))
    end
end

@inline Base.zero(::Mask{W,U}) where {W,U} = Mask{W}(zero(U))

@inline data(m::Mask) = m.u
@generated function Base.:(&)(m1::Mask{W,U}, m2::Mask{W,U}) where {W,U}
    binary_mask_op(W, U, "and")
end
@generated function Base.:(|)(m1::Mask{W,U}, m2::Mask{W,U}) where {W,U}
    binary_mask_op(W, U, "or")
end
@generated function Base.:(⊻)(m1::Mask{W,U}, m2::Mask{W,U}) where {W,U}
    binary_mask_op(W, U, "xor")
end
@generated function Base.:(==)(m1::Mask{W,U}, m2::Mask{W,U}) where {W,U}
    binary_mask_op(W, U, "icmp eq")
end
@generated function Base.:(!=)(m1::Mask{W,U}, m2::Mask{W,U}) where {W,U}
    binary_mask_op(W, U, "icmp ne")
end

@generated function splitint(i::S, ::Type{T}) where {S <: Base.BitInteger, T <: Union{Bool,Base.BitInteger}}
    sizeof_S = sizeof(S)
    sizeof_T = sizeof(T)
    if sizeof_T > sizeof_S
        return :(i % T)
    elseif sizeof_T == sizeof_S
        return :i
    end
    W, r = divrem(sizeof_S, sizeof_T)
    @assert iszero(r)
    vtyp = "<$W x i$(8sizeof_T)>"
    instrs = """
        %split = bitcast i$(8sizeof_S) %0 to $vtyp
        ret $vtyp %split
    """
    quote
        $(Expr(:meta,:inline))
        Vec(llvmcall($instrs, _Vec{$W,$T}, Tuple{$S}, i))
    end
end
@generated function fuseint(v::Vec{W,I}) where {W, I <: Union{Bool,Base.BitInteger}}
    @assert ispow2(W)
    bytes = W * sizeof(I)
    bits = 8bytes
    @assert bytes ≤ 16
    T = (I <: Signed) ? Symbol(:Int, bits) : Symbol(:UInt, bits)
    vtyp = "<$W x i$(8sizeof(I))>"
    styp = "i$(bits)"
    instrs = """
        %fused = bitcast $vtyp %0 to $styp
        ret $styp %fused
    """
    quote
        $(Expr(:meta,:inline))
        llvmcall($instrs, $T, Tuple{_Vec{$W,$I}}, data(v))
    end
end

function vadd_expr(W,U)
    instrs = String[]
    truncate_mask!(instrs, '0', W, 0)
    truncate_mask!(instrs, '1', W, 1)
    push!(instrs, """%uv.0 = zext <$W x i1> %mask.0 to <$W x i8>
    %uv.1 = zext <$W x i1> %mask.1 to <$W x i8>
    %res = add <$W x i8> %uv.0, %uv.1
    ret <$W x i8> %res""")
    Expr(:block, Expr(:meta, :inline), :(Vec(llvmcall($(join(instrs, "\n")), _Vec{$W,UInt8}, Tuple{$U, $U}, m1.u, m2.u))))
end
@generated Base.:(+)(m1::Mask{W,U}, m2::Mask{W,U}) where {W,U} = vadd_expr(W,U)

@inline Base.:(&)(m::Mask{W}, b::Bool) where {W} = Mask{W}(b ? m.u : zero(m.u))
@inline Base.:(&)(b::Bool, m::Mask{W}) where {W} = Mask{W}(b ? m.u : zero(m.u))

@inline Base.:(|)(m::Mask{W,U}, b::Bool) where {W,U} = b ? max_mask(Mask{W,U}) : m
@inline Base.:(|)(b::Bool, m::Mask{W,U}) where {W,U} = b ? max_mask(Mask{W,U}) : m

@inline Base.:(⊻)(m::Mask{W}, b::Bool) where {W} = Mask{W}(b ? ~m.u : m.u)
@inline Base.:(⊻)(b::Bool, m::Mask{W}) where {W} = Mask{W}(b ? ~m.u : m.u)

@inline Base.:(<<)(m::Mask{W}, i::IntegerTypesHW) where {W} = Mask{W}(shl(m.u, i))
@inline Base.:(>>)(m::Mask{W}, i::IntegerTypesHW) where {W} = Mask{W}(shr(m.u, i))
@inline Base.:(>>>)(m::Mask{W}, i::IntegerTypesHW) where {W} = Mask{W}(shr(m.u, i))

for (U,W) in [(UInt8,8), (UInt16,16), (UInt32,32), (UInt64,64)]
    @eval @inline vany(m::Mask{$W,$U}) = m.u != $(zero(U))
    @eval @inline vall(m::Mask{$W,$U}) = m.u == $(typemax(U))
end
@inline vany(m::Mask{W}) where {W} = (m.u & max_mask(Val{W}()).u) !== zero(m.u)
@inline vall(m::Mask{W}) where {W} = (m.u & max_mask(Val{W}()).u) === (max_mask(Val{W}()).u)

@generated function Base.:(!)(m::Mask{W,U}) where {W,U}
    mtyp_input = "i$(8sizeof(U))"
    mtyp_trunc = "i$(W)"
    instrs = String[]
    truncate_mask!(instrs, '0', W, 0)
    mask = llvmconst(W, "i1 true")
    push!(instrs, "%resvec.0 = xor <$W x i1> %mask.0, $mask")
    zext_mask!(instrs, "resvec.0", W, 1)
    push!(instrs, "ret $mtyp_input %res.1")
    quote
        $(Expr(:meta,:inline))
        Mask{$W}(llvmcall($(join(instrs,"\n")), $U, Tuple{$U}, m.u))
    end
end
@inline Base.:(~)(m::Mask) = !m

@inline Base.count_ones(m::Mask) = count_ones(m.u)
@inline Base.:(+)(m::Mask, i::Integer) = i + count_ones(m)
@inline Base.:(+)(i::Integer, m::Mask) = i + count_ones(m)

function mask_type_symbol(W)
    if W <= 8
        return :UInt8
    elseif W <= 16
        return :UInt16
    elseif W <= 32
        return :UInt32
    elseif W <= 64
        return :UInt64
    else#if W <= 128
        return :UInt128
    end
end
function mask_type(W)
    if W <= 8
        return UInt8
    elseif W <= 16
        return UInt16
    elseif W <= 32
        return UInt32
    elseif W <= 64
        return UInt64
    else#if W <= 128
        return UInt128
    end
end
mask_type(::Union{Val{1},StaticInt{1}}) = UInt8#Bool
mask_type(::Union{Val{2},StaticInt{2}}) = UInt8
mask_type(::Union{Val{4},StaticInt{4}}) = UInt8
mask_type(::Union{Val{8},StaticInt{8}}) = UInt8
mask_type(::Union{Val{16},StaticInt{16}}) = UInt16
mask_type(::Union{Val{32},StaticInt{32}}) = UInt32
mask_type(::Union{Val{64},StaticInt{64}}) = UInt64

@generated function mask_type(::Type{T}, ::Union{Val{P},StaticInt{P}}) where {T,P}
    mask_type_symbol(pick_vector_width(P, T))
end
@generated function mask_type(::Type{T}) where {T}
    W = max(1, register_size(T) >>> intlog2(T))
    mask_type_symbol(W)
    # mask_type_symbol(pick_vector_width(T))
end
@generated function Base.zero(::Type{<:Mask{W}}) where {W}
    Expr(:block, Expr(:meta, :inline), Expr(:call, Expr(:curly, :Mask, W), Expr(:call, :zero, mask_type_symbol(W))))
end

@generated function max_mask(::Union{Val{W},StaticInt{W}}) where {W}
    U = mask_type(W)
    Mask{W,U}(one(U)<<W - one(U))
end
@inline max_mask(::Type{T}) where {T} = max_mask(pick_vector_width_val(T))
@generated max_mask(::Type{Mask{W,U}}) where {W,U} = Mask{W,U}(one(U)<<W - one(U))

@generated function valrem(::Union{Val{W},StaticInt{W}}, l) where {W}
    ex = ispow2(W) ? :(l & $(W - 1)) : :(l % $W)
    Expr(:block, Expr(:meta, :inline), ex)
end
@generated function mask(::Union{Val{W},StaticInt{W}}, l::Integer) where {W}
    M = mask_type(W)
    quote
        $(Expr(:meta,:inline))
        rem = valrem(Val{$W}(), vsub((l % $M), one($M)))
        Mask{$W,$M}($(typemax(M)) >>> ($(M(8sizeof(M))-1) - rem))
    end
end
@generated mask(::Union{Val{W},StaticInt{W}}, ::StaticInt{L}) where {W, L} = mask(Val(W), L)
@inline mask(::Type{T}, l::Integer) where {T} = mask(pick_vector_width_val(T), l)

# @generated function masktable(::Union{Val{W},StaticInt{W}}, rem::Integer) where {W}
#     masks = Expr(:tuple)
#     for w ∈ 0:W-1
#         push!(masks.args, data(mask(Val(W), w == 0 ? W : w)))
#     end
#     Expr(
#         :block,
#         Expr(:meta,:inline),
#         Expr(:call, Expr(:curly, :Mask, W), Expr(
#             :macrocall, Symbol("@inbounds"), LineNumberNode(@__LINE__, Symbol(@__FILE__)),
#             Expr(:call, :getindex, masks, Expr(:call, :+, 1, Expr(:call, :valrem, Expr(:call, Expr(:curly, W)), :rem)))
#         ))
#     )
# end

@inline tomask(m::Unsigned) = m
@inline tomask(m::Mask) = m
@generated function tomask(v::Vec{W,Bool}) where {W}
    usize = W > 8 ? nextpow2(W) : 8
    utyp = "i$(usize)"
    U = mask_type_symbol(W)
    instrs = String[]
    push!(instrs, "%bitvec = trunc <$W x i8> %0 to <$W x i1>")
    zext_mask!(instrs, "bitvec", W, 0)
    push!(instrs, "ret i$(usize) %res.0")
    quote
        $(Expr(:meta, :inline))
        Mask{$W}(llvmcall($(join(instrs, "\n")), $U, Tuple{_Vec{$W,Bool}}, data(v)))
    end
end
@inline tomask(v::AbstractSIMDVector{<:Any,Bool}) = tomask(data(v))

@generated function Base.:(%)(m::Mask{W,U}, ::Type{I}) where {W,U,I<:Integer}
    bits = 8sizeof(I)
    instrs = String[]
    truncate_mask!(instrs, '0', W, 0)
    push!(instrs, "%res = zext <$W x i1> %mask.0 to <$W x i$(bits)>\nret <$W x i$(bits)> %res")
    quote
        $(Expr(:meta,:inline))
        Vec(llvmcall($(join(instrs, "\n")), _Vec{$W,$I}, Tuple{$U}, data(m)))
    end
end
Vec(m::Mask{W}) where {W} = m % int_type(Val{W}())

# @inline getindexzerobased(m::Mask, i) = (m.u >>> i) % Bool
# @inline function extractelement(m::Mask{W}, i::Integer) where {W}
#     @boundscheck i > W && throw(BoundsError(m, i))
#     getindexzerobased(m, i)
# end
@generated function extractelement(v::Mask{W,U}, i::I) where {W,U,I}
    instrs = String[]
    truncate_mask!(instrs, '0', W, 0)
    push!(instrs, "%res1 = extractelement <$W x i1> %mask.0, i$(8sizeof(I)) %1")
    push!(instrs, "%res8 = zext i1 %res1 to i8\nret i8 %res8")
    instrs_string = join(instrs, "\n")
    call = :(llvmcall($instrs_string, Bool, Tuple{$U,$I}, data(v), i))
    Expr(:block, Expr(:meta, :inline), call)
end
@generated function insertelement(v::Mask{W,U}, x::T, i::I) where {W, T, U, I <: Union{Bool,IntegerTypesHW}}
    mtyp_input = "i$(max(8,W))"
    instrs = String["%bit = trunc i$(8sizeof(T)) %1 to i1"]
    truncate_mask!(instrs, '0', W, 0)
    push!(instrs, "%bitvec = insertelement <$W x i1> %mask.0, i1 %bit, i$(8sizeof(I)) %2")
    zext_mask!(instrs, "bitvec", W, 1)
    push!(instrs, "ret $(mtyp_input) %res.1")
    instrs_string = join(instrs, "\n")
    call = :(Mask{$W}(llvmcall($instrs_string, $U, Tuple{$U,$T,$I}, data(v), x, i)))
    Expr(:block, Expr(:meta, :inline), call)
end

# @generated function Base.isodd(i::MM{W,1}) where {W}
#     U = mask_type(W)
#     evenfirst = 0xaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaa % U
#     # Expr(:block, Expr(:meta, :inline), :(isodd(i.i) ? Mask{$W}($oddfirst) : Mask{$W}($evenfirst)))
#     Expr(:block, Expr(:meta, :inline), :(Mask{$W}($evenfirst >> (i.i & 0x03))))
# end
# @generated function Base.iseven(i::MM{W,1}) where {W}
#     U = mask_type(W)
#     oddfirst = 0x55555555555555555555555555555555 % U
#     # evenfirst = oddfirst << 1
#     # Expr(:block, Expr(:meta, :inline), :(isodd(i.i) ? Mask{$W}($evenfirst) : Mask{$W}($oddfirst)))
#     Expr(:block, Expr(:meta, :inline), :(Mask{$W}($oddfirst >> (i.i & 0x03))))
# end
@inline Base.isodd(i::MM{W,1}) where {W} = Mask{W}((0xaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaa % mask_type(Val{W}())) >>> (i.i & 0x03))
@inline Base.iseven(i::MM{W,1}) where {W} = Mask{W}((0x55555555555555555555555555555555 % mask_type(Val{W}())) >>> (i.i & 0x03))

function cmp_quote(W, cond, vtyp, T1, T2 = T1)
    instrs = String["%m = $cond $vtyp %0, %1"]
    zext_mask!(instrs, 'm', W, '0')
    push!(instrs, "ret i$(max(8,W)) %res.0")
    U = mask_type(W);
    quote
        $(Expr(:meta,:inline))
        Mask{$W}(llvmcall($(join(instrs, "\n")), $U, Tuple{_Vec{$W,$T1},_Vec{$W,$T2}}, data(v1), data(v2)))
    end
end
function icmp_quote(W, cond, bytes, T1, T2 = T1)
    vtyp = vtype(W, "i$(8bytes)");
    cmp_quote(W, "icmp " * cond, vtyp, T1, T2)
end
function fcmp_quote(W, cond, T)
    vtyp = vtype(W, T === Float32 ? "float" : "double");
    cmp_quote(W, "fcmp nsz arcp contract reassoc " * cond, vtyp, T)
end
# @generated function compare(::Val{cond}, v1::Vec{W,I}, v2::Vec{W,I}) where {cond, W, I}
    # cmp_quote(W, cond, sizeof(I), I)
# end
# for (f,cond) ∈ [(:(==), :eq), (:(!=), :ne), (:(>), :ugt), (:(≥), :uge), (:(<), :ult), (:(≤), :ule)]
for (f,cond) ∈ [(:(==), "eq"), (:(!=), "ne")]
    @eval @generated function Base.$f(v1::Vec{W,T1}, v2::Vec{W,T2}) where {W,T1<:Integer,T2<:Integer}
        @assert sizeof(T1) == sizeof(T2)
        icmp_quote(W, $cond, sizeof(T1), T1, T2)
    end
end
for (f,cond) ∈ [(:(>), "ugt"), (:(≥), "uge"), (:(<), "ult"), (:(≤), "ule")]
    @eval @generated function Base.$f(v1::Vec{W,T}, v2::Vec{W,T}) where {W,T<:Unsigned}
        icmp_quote(W, $cond, sizeof(T), T)
    end
    # bytes = 1
    # while bytes ≤ 8
    #     T = Symbol(:UInt, 8bytes)
    #     W = 2
    #     while W * bytes ≤ REGISTER_SIZE
    #         @eval @inline Base.$f(v1::Vec{$W,$T}, v2::Vec{$W,$T}) = $(icmp_quote(W, cond, bytes, T))
    #         W += W
    #     end
    #     bytes += bytes
    # end
    # @eval @inline Base.$f(v1::Vec{W,I}, v2::Vec{W,I}) where {W, I <: Signed} = compare(Val{$(QuoteNode(cond))}(), v1, v2)
end
for (f,cond) ∈ [(:(>), "sgt"), (:(≥), "sge"), (:(<), "slt"), (:(≤), "sle")]
    @eval @generated function Base.$f(v1::Vec{W,T}, v2::Vec{W,T}) where {W,T<:Signed}
        icmp_quote(W, $cond, sizeof(T), T)
    end
    # bytes = 1
    # while bytes ≤ 8
    #     T = Symbol(:Int, 8bytes)
    #     W = 2
    #     while W * bytes ≤ REGISTER_SIZE
    #         @eval @inline Base.$f(v1::Vec{$W,$T}, v2::Vec{$W,$T}) = $(cmp_quote(W, cond, bytes, T))
    #         W += W
    #     end
    #     bytes += bytes
    # end
end

# for (f,cond) ∈ [(:(==), "oeq"), (:(>), "ogt"), (:(≥), "oge"), (:(<), "olt"), (:(≤), "ole"), (:(≠), "one")]
for (f,cond) ∈ [(:(==), "ueq"), (:(>), "ugt"), (:(≥), "uge"), (:(<), "ult"), (:(≤), "ule"), (:(≠), "une")]
    @eval @generated function Base.$f(v1::Vec{W,T}, v2::Vec{W,T}) where {W, T <: Union{Float32,Float64}}
        fcmp_quote(W, $cond, T)
    end
    # for (T,bytes) ∈ [(:Float32,4), (:Float64,8)]
    #     W = 2
    #     while W * bytes ≤ REGISTER_SIZE
    #         @eval @inline Base.$f(v1::Vec{$W,$T}, v2::Vec{$W,$T}) = $(fcmp_quote(W, cond, T))
    #         W += W
    #     end
    # end
end

# import IfElse: ifelse
@generated function ifelse(m::Mask{W,U}, v1::Vec{W,T}, v2::Vec{W,T}) where {W,U<:Unsigned,T}
    typ = LLVM_TYPES[T]
    vtyp = vtype(W, typ)
    selty = vtype(W, "i1")
    f = "select"
    if Base.libllvm_version ≥ v"9" && ((T === Float32) || (T === Float64))
        f *= " nsz arcp contract reassoc"
    end
    instrs = String[]
    truncate_mask!(instrs, '0', W, 0)
    push!(instrs, "%res = $f $selty %mask.0, $vtyp %1, $vtyp %2\nret $vtyp %res")
    quote
        $(Expr(:meta,:inline))
        Vec(llvmcall($(join(instrs,"\n")), _Vec{$W,$T}, Tuple{$U,_Vec{$W,$T},_Vec{$W,$T}}, data(m), data(v1), data(v2)))
    end
end
# @inline ifelse(m::Mask, v::Vec, s) = ((x,y) = promote(v,s); ifelse(m,x,y))
# @inline ifelse(m::Mask, s, v::Vec) = ((x,y) = promote(s,v); ifelse(m,x,y))
@inline ifelse(m::Mask{W}, s1::T, s2::T) where {W,T<:NativeTypes} = ifelse(m, Vec{W,T}(s1), Vec{W,T}(s2))
@inline ifelse(m::Mask{W}, s1, s2) where {W} = ((x1,x2) = promote(s1,s2); ifelse(m, x1, x2))

@inline Base.isnan(v::AbstractSIMD) = v != v

@inline Base.isfinite(x::AbstractSIMD) = iszero(x - x)

@inline Base.flipsign(x::AbstractSIMD, y::AbstractSIMD) = ifelse(y > zero(y), x, -x)
for T ∈ [:Float32, :Float64]
    @eval begin
        @inline Base.flipsign(x::AbstractSIMD, y::$T) = ifelse(y > zero(y), x, -x)
        @inline Base.flipsign(x::$T, y::AbstractSIMD) = ifelse(y > zero(y), x, -x)
    end
end
@inline Base.flipsign(x::AbstractSIMD, y::Real) = ifelse(y > zero(y), x, -x)
@inline Base.flipsign(x::Real, y::AbstractSIMD) = ifelse(y > zero(y), x, -x)
@inline Base.flipsign(x::Signed, y::AbstractSIMD) = ifelse(y > zero(y), x, -x)
@inline Base.isodd(x::AbstractSIMD{W,T}) where {W,T<:Integer} = (x & one(T)) != zero(T)

@generated function ifelse(m::Vec{W,Bool}, v1::Vec{W,T}, v2::Vec{W,T}) where {W,T}
    typ = LLVM_TYPES[T]
    vtyp = vtype(W, typ)
    selty = vtype(W, "i1")
    f = "select"
    if Base.libllvm_version ≥ v"9" && ((T === Float32) || (T === Float64))
        f *= " nsz arcp contract reassoc"
    end
    instrs = String["%mask.0 = trunc <$W x i8> %0 to <$W x i1>"]
    # truncate_mask!(instrs, '0', W, 0)
    push!(instrs, "%res = $f $selty %mask.0, $vtyp %1, $vtyp %2\nret $vtyp %res")
    quote
        $(Expr(:meta,:inline))
        Vec(llvmcall($(join(instrs,"\n")), _Vec{$W,$T}, Tuple{_Vec{$W,Bool},_Vec{$W,$T},_Vec{$W,$T}}, data(m), data(v1), data(v2)))
    end
end
@inline ifelse(b::Bool, s::NativeTypes, v::V) where {V <: AbstractSIMD} = ifelse(b, convert(V, s), v)
@inline ifelse(b::Bool, v::V, s::NativeTypes) where {V <: AbstractSIMD} = ifelse(b, v, convert(V, s))

@generated function Base.convert(::Type{Bit}, v::Vec{W,Bool}) where {W,Bool}
    instrs = String[]
    push!(instrs, "%m = trunc <$W x i8> %0 to <$W x i1>")
    zext_mask!(instrs, 'm', W, '0')
    push!(instrs, "ret i$(max(8,W)) %res.0")
    U = mask_type(W);
    quote
        $(Expr(:meta,:inline))
        Mask{$W}(llvmcall($(join(instrs, "\n")), $U, Tuple{_Vec{$W,Bool}}, data(v)))
    end
end
@inline Base.convert(::Type{Vec{W,Bit}}, v::Vec{W,Bool}) where {W,Bool} = convert(Bit, v)


text:
# This file was generated, do not modify it. # hide
r_depth = range(dtr_model, :max_depth, lower=2, upper=20)

tm = TunedModel(
    model=dtr_model,
    ranges=[r_depth],
    tuning=Grid(resolution=10),
    resampling=CV(nfolds=5, rng=1254),
    measure=rms
)

mtm = machine(tm, X, y)

fit!(mtm, rows=train)

ypred = MLJ.predict(mtm, rows=test)
round(rms(ypred, y[test]), sigdigits=3)

max_issues_repo_name: alan-turing-institute/MLJTutorials

max_forks_repo_name: alan-turing-institute/MLJTutorials


text:
export where

immutable VectorAlgos <: AFImpl
	where::Ptr{Void}

	function VectorAlgos(ptr)
		new(
			Libdl.dlsym(ptr, :af_where)
		)
	end
end

@afCall_Arr_Arr(where, vectorAlgos, where)

max_stars_repo_name: unbornchikken/julia-ml-proto

max_issues_repo_name: unbornchikken/julia-ml-proto

max_forks_repo_name: unbornchikken/julia-ml-proto


text:
function allunitsharvest!(game)
	#groundlevel=false
	#if unit.pl=[2]
	#	groundlevel=true
	#end
	#bools=(unit.passover,unit.passoverself,unit.inclusive)
	#influencemap=allinfluence(game,unit.groundlevel,bools)
	allunitslive!(game)
	points=[0.0,0,0,0,0]
	for (loc,unit) in game.map
		if unit!=0
			points.+=unit.harvest!(game,unit)
#			println(points)
		end
	end
	#game.points.+=points
	game.season+=1
	push!(game.sequence,:harvest)
	GAccessor.text(game.g[1,2],pointslabel(game))
	return points
end

function lifluence(game,unit::Unit)
	lifemap=unitslive(game,unit.color)
	connectedunits=Unit[unit]
	cwhite=Unit[]
	group=[unit.loc]
	if unit.color==(1,1,1)
		push!(cwhite,unit)
		stuff=Dict()
		#stuff[:loced]=[unit.loc]
		reach=makegrid(7)
		for loc in reach
			u=game.map[loc]
			if isa(u,Unit) && u.color!=(1,1,1) && !haskey(stuff,u.color)	#!in(loc,stuff[:loced])
				lmu=lifluence(game,u)
				if in(unit,lmu[3])
					stuff[u.color]=lmu
				end
			end
		end
		for (col,lm) in stuff
			for l in lm[1]
				if !in(l,group)
					push!(group,l)
				end
			end
			for cu in lm[2]
				if !in(cu,connectedunits)
					push!(connectedunits,cu)
				end
			end
			for cw in lm[3]
				if !in(cw,cwhite)
					push!(cwhite,cw)
				end
			end
		end
		return (group,connectedunits,cwhite)
	end
	temp=[unit.loc]
	while !isempty(temp)
		temp2=Tuple[]
		for t in temp
			for h in adjacent(t)
#				lmt=lifemap[t]
#				lmh=lifemap[h]
#				canspread=true
				if !in(h,group) && !in(h,temp) && !in(h,temp2) && in(h,game.board.grid) && sum(lifemap[h])>0
					push!(temp2,h)
					u=game.map[h]
					if isa(u,Unit)
						push!(connectedunits,u)
						if u.color==(1,1,1)
							push!(cwhite,u)
						end
					end
				end
			end
		end
		for t2 in temp2
			push!(group,t2)
		end
		temp=temp2
	end
	return (group,connectedunits,cwhite)
end
function connectedunits(game,unit)
	(lif,cu,cw)=lifluence(game,unit)
	cellconnected=getcellgroup(game,unit)
	return (cw,cellconnected,cu) #connected white units, connected by cells, connected by lifluence
end

function influence(game,hex,groundlevel=false,passover=false,passoverself=true,inclusive=true)
	unit=game.map[hex]
	white=(1,1,1)
	group=Dict(hex=>6.0)
	temp=Dict(hex=>6.0)
	for rad in 1:unit.ir
		temp2=Dict()
		for t in temp
			for h in adjacent(t[1],1,groundlevel)
				if !in(h,keys(group)) && !in(h,keys(temp)) && !in(h,keys(temp2)) && in(h,keys(game.map))
					inf=1/rad
					if game.map[h]==0 || passover
						temp2[h]=inf
					elseif passoverself && (game.map[h].color==unit.color || game.map[h].color==white)
						temp2[h]=inf
					end
					if inclusive && game.map[h]!=0 && !in(h,keys(temp2))
						group[h]=inf
					end
				end
			end
		end
		for (h2,i2) in temp2
			group[h2]=i2
		end
		temp=temp2
	end
	return group
end
function allinfluence(game,groundlevel=false,bools=(true,false,true,true))
	influencemap=Dict()
	for (loc,player) in game.map
		if !groundlevel || loc[3]==2
			influencemap[loc]=[0.0,0,0]
		end
	end
	for (loc,unit) in game.map
		if unit!=0
			col=unit.color
			infl=influence(game,loc,bools...)
			for inf in infl
				influencemap[inf[1]].+=inf[2].*col
			end
		end
	end
	return influencemap
end

function peekharvest(game,groundlevel=false,bools=(true,false,true,true))
	influencemap=allinfluence(game,groundlevel,bools)
	brgbw=[0.0,0,0,0,0]
	for (iloc,inf) in influencemap
		ninf=numcolors(inf)
		if ninf==3
			brgbw[5]+=min(inf...)
		end
		if ninf==1
			brgbw[1]+=min(sum(inf),1)
		else
			for c in 1:3
				brgbw[c+1]+=min(inf[c],inf[c%3+1])
			end
		end
	end
	for c in 2:4
		brgbw[c]-=brgbw[5]
	end
	return brgbw
end

function getgroup_dep(game,hex)
	player=game.map[hex]
	white=(1,1,1)
	if player==0
		return []
	end
	group=Tuple[hex]
	temp=[hex]
	while !isempty(temp)
		temp2=Tuple[]
		for t in temp
			for h in adjacent(t)
				if !in(h,group) && !in(h,temp) && !in(h,temp2) && in(h,keys(game.map)) && (game.map[h]==player || game.map[h]==white)
					push!(temp2,h)
				end
			end
		end
		for t2 in temp2
			push!(group,t2)
		end
		temp=temp2
	end
	return group
end

function unitcost(game,loc,unitparams)
	distance=nearestwhite(game,loc,loc[3]==2)
	cost=distance*unitparams[2]
	cost*=sum(unitparams[1]) #only if sum>1?
	return cost
end
function unitcost(game,unit::Unit)
	return unit.costfun(game)
end
function subtractcost(game,cost,color)
	rgb=cost.*color
	srgb=sum(rgb)
	if rgb[1]<=game.points[2] && rgb[2]<=game.points[3] && rgb[3]<=game.points[4]
		game.points[2]-=rgb[1]
		game.points[3]-=rgb[2]
		game.points[4]-=rgb[3]
	elseif srgb/2<=game.points[5] && srgb/2<=game.points[1]
		game.points[1]-=srgb/2
		game.points[5]-=srgb/2
	else
		return false
	end
	return true
end
function updategroups_dep!(game::Game) #sometimes doesn't find most recent unit..? Or is there an issue when deleting units? Maybe fixed. Now it sometimes adds the most recent unit twice...
	groups=Group[]
	for spawn in game.spawns
		push!(groups,getgroup(game,spawn))
	end
	unique=Group[groups[1]]
	for group in groups
		for uniq in unique
			if samegroup(group,uniq)
				break
			end
			push!(unique,group)
		end
	end
	game.groups=unique
end

function getgroup_dep(game,unit::Unit,color=-1,connectedunits=Unit[]) #why don't white units get added to spawns? Maybe they shouldn't, bug in our favor. They should be available as partial spawns, divide white into 3 spawn units. Fixed some stuff, still good bugs? Well, this should be rewritten rather than bugged down, first a initgroup that gets the body.
	if color==-1
		color=unit.color
	end
	lifemap=unitslive(game,unit.color)
	if !in(unit,connectedunits)
		push!(connectedunits,unit)
	end
	cwhite=Unit[] #connected spawns
	if color==(1,1,1) #this needs to be rewritten to allow colored spawns, sorta works now but... nonbody groups automerge with white group
		push!(cwhite,unit)
		stuff=Dict()
		reach=makegrid(7) #should be minimum twice the maximum ir of all units +1, should be a better way
		for lo in reach
			loc=lo.+(unit.loc.-(0,0,2))
			if !in(loc,keys(game.map))
				continue
			end
			u=game.map[loc]
			if isa(u,Unit) && distance(u.loc,unit.loc)<(u.ir+unit.ir+2) && !haskey(stuff,u.color) && !in(u,connectedunits) && u.color!=(1,1,1)
				subgroup=getgroup(game,u,u.color,connectedunits)
				if in(unit,subgroup.spawns)
					stuff[u.color]=subgroup
				end
			end
		end
		for (col,sg) in stuff
			for cu in sg.units
				if !in(cu,connectedunits)
					push!(connectedunits,cu)
				end
			end
			for cw in sg.spawns
				if !in(cw,cwhite)
					push!(cwhite,cw)
				end
			end
		end
		body=getcellgroup(game,cwhite[1])
		return newgroup(cwhite,body,connectedunits)
	end
	ulocs=[unit.loc]
	temp=[unit.loc]
	while !isempty(temp)
		temp2=Tuple[]
		for t in temp
			for h in adjacent(t)
				if !in(h,ulocs) && !in(h,temp) && !in(h,temp2) && in(h,game.board.grid) && sum(lifemap[h])>0
					push!(temp2,h)
					u=game.map[h]
					if isa(u,Unit) && !in(u,connectedunits)
						push!(connectedunits,u)
						if u.color==(1,1,1)
							push!(cwhite,u)
							for adjwu in adjacent(u.loc)
								adju=game.map[adjwu]
								if isa(adju,Unit) && adju.color==color && !in(adju,connectedunits)
									nsg=getgroup(game,adju,color,connectedunits)
									for cu in nsg.units
										if !in(cu,connectedunits)
											push!(connectedunits,cu)
										end
									end
									for cw in nsg.spawns
										if !in(cw,cwhite)
											push!(cwhite,cw)
										end
									end
								end
							end
						end
					end
				end
			end
		end
		for t2 in temp2
			push!(ulocs,t2)
		end
		temp=temp2
	end
	body=getcellgroup(game,connectedunits[1]) #tempfix
	return newgroup(cwhite,body,connectedunits)
end

function allunitslive!(game)
	lifemap=Dict()
	for loc in game.board.grid
		lifemap[loc]=[0.0,0,0]
	end
	for (loc,unit) in game.map
		if isa(unit,Unit)
			unit.live!(game,unit,lifemap)
		end
	end
	game.lifemap=lifemap
	return lifemap
end
#lifemap=allunitslive! #maybe not store lifemap in game since it gets modified in so many places. Rewrite everything to work with a temp lifemap. Or rather have a solid lifemap and a temp harvest ledger

function printpoints(game)
	harv=round.(peekharvest(game),1,10)
	println("Black: ",harv[1]," Red: ",harv[2]," Green: ",harv[3]," Blue: ",harv[4]," White: ",harv[5])#," Total: ",sum(harv))
end
function nearestwhite(game,hex,layer=false) #rewrite to check for spawn? Or spawning is always done from groups? This isn't even needed now that costs are gone
	white=(1,1,1)
	group=[hex]
	temp=[hex]
	for rad in 1:1000000
		temp2=[]
		for t in temp
			for h in adjacent(t,1,layer)
				if in(h,keys(game.map))
					unit=game.map[h]
					if unit!=0 && unit.color==white
						return rad
					end
					if !in(h,group) && !in(h,temp) && !in(h,temp2)
						push!(temp2,h)
					end
				end
			end
		end
		for h2 in temp2
			push!(group,h2)
		end
		temp=temp2
	end
	return Inf
end
function getexpandcost(shells::Integer=6,initlocs=[(6,6,2)],basecost=50) #deprecated
	patch=makegrid(shells,initlocs)
	cost=length(patch)+basecost*length(initlocs)
	return cost
end
function harvest!(game,group::Group)
	if !group.harvested
		harv=checkharvest(game,group) #some stuff needs to be rewritten to avoid the deepcopy above
		group.points+=harv
		group.harvested=true
		for unit in group.units
			unit.harvested=true
		end
	end
end
function allgroupsharvest!(game)
	for group in game.groups
		harvest!(game,group)
	end
end

function collectharvest!(game::Game)
	for group in game.groups
		game.points+=group.points
		group.points-=group.points
	end
end
function newseason!(game)
	for group in game.groups
		group.harvested=false
		for unit in group.units
			unit.harvested=false
		end
	end
	game.season+=1
end
function harvest!(game::Game)
	allunitslive!(game) #maybe call only when new units are placed. But then placing units takes longer...
	updategroups!(game)
	allgroupsharvest!(game)
	collectharvest!(game)
	newseason!(game)
	push!(game.sequence,:harvest)
	GAccessor.text(game.g[1,2],pointslabel(game))
	return game.points
end


text:
struct DataOneV1 <: DataRepo
end

# The only APIs known to be supported by DataOne Version 1 is DataDryad
# As and when we come to know of other APIs supported by DataOne, more abstraction layers will be introduced.
# Named as DataOneV1 to make new contributors to the package aware that DataOneV1 exists.
base_url(::DataOneV1) = "https://datadryad.org/mn/object/http://dx.doi.org/"

description(repo::DataOneV1, mainpage) = text_only(first(eachmatch(sel"dcterms\:description", mainpage.root)))

author(::DataOneV1, mainpage) = [text_only(i) for i in eachmatch(sel"dcterms\:creator", mainpage.root)]

license(::DataOneV1, mainpage) = "http://creativecommons.org/publicdomain/zero/1.0/"

published_date(::DataOneV1, mainpage) = Dates.DateTime(chop(text_only(first(eachmatch(sel"dcterms\:dateSubmitted", mainpage.root)))))

paper_cite(::DataOneV1, mainpage) = citation_text(text_only(first(eachmatch(sel"dcterms\:references", mainpage.root))))

dataset_cite(::DataOneV1, mainpage) = citation_text(text_only(first(eachmatch(sel"dcterms\:identifier", mainpage.root))))

function get_urls(repo::DataOneV1, page)
    urls = []
    links = eachmatch(sel"dcterms\:hasPart", page.root)
    for link in links
        downloadlinks = eachmatch(sel".image-link", getpage(text_only(link)).root)
        for downloadlink in downloadlinks
            push!(urls, downloadlink.attributes["href"])
        end
    end
    urls
end

function get_checksums(repo::DataOneV1, page)
    checksums = []
    links = eachmatch(sel"dcterms\:hasPart", page.root)
    for link in links
        spans = eachmatch(sel".file-list span", getpage(text_only(link)).root)
        for span in spans
            if occursin(r"^[a-f0-9]{32}$", text_only(span))
                push!(checksums, (:md5, text_only(span)))
            end
        end
    end
    if length(checksums) > 0
        @info("The generated registration block uses md5 hash, " *
            "the MD5.jl package must be loaded to run the registration")
    end
    checksums
end

function data_fullname(::DataOneV1, mainpage)
    text_only(first(eachmatch(sel"dcterms\:title", mainpage.root)))
end

function website(::DataOneV1, mainpage_url, mainpage)
    replace(mainpage_url, "https://datadryad.org/mn/object/" => "")
end


text:
using JuliaFormatter
format(".")


text:
using TropicalCycloneRecords

atlantic_url ="https://www.nhc.noaa.gov/data/hurdat/hurdat2-1851-2019-042820.txt"
data = Dict{String,Dict{T,N} where T where N}()
res = Tropical.read_hurdat(atlantic_url,true)
Tropical.hurdat_to_dictionary(res,data)


text:
#
# This file is part of the Actors.jl Julia package,
# MIT license, part of https://github.com/JuliaActors
#

using Actors
import Actors: spawn

function fib(D::Dict{Int,BigInt}, n::Int)
    get!(D, n) do
        n == 0 && return big(0)
        n == 1 && return big(1)
        return fib(D, n-1) + fib(D, n-2)
    end
end

myfib = spawn(Bhv(fib, Dict{Int,BigInt}()))

request(myfib, 1000)


text:
# using DataFrames
# showln(x) = (show(x); println())
#
# # A DataFrame is an in-memory database
# df = DataFrame(A = [1, 2], B = [ℯ, π], C = ["xx", "xy"])
# showln(df)

cd(dirname(@__FILE__()))

using Plots
using PlotThemes
using Measures
using LinearAlgebra
using SparseArrays

#=

Part I: A Simple Implementation of Newton's Method

=#

################################################################################
# Defining Constructors (Types)
################################################################################

#=
We wish to solve F(x) == 0. For this we need to
- Evaluate F(x)
- Approximate or Calcualte F'(x)
- Report || F(x) || : For some problems, ||⋅|| is not the Euclidan norm

We thus define a simple constructor (object) with three fields:
=#
mutable struct Residual
    evalF::Function     # F(x)
    jac::Function       # F'(x)
    Fnrm::Function      # || F(x) ||
end

# We make Residual a callable type that acts as a function but encapsulates
# data or other important functions.
function (F::Residual)(x::Array{Float64,1})
    return F.evalF(x)
end

################################################################################
# Basic Newton's Method
################################################################################

# This is the basic structure of a classic version of Newton's method for
# nonlinear equations. We will slightly adapt it blow to fit the applications.
function newton_nd(F::Residual,x_0::Array{Float64,1},
                   atol::Float64,rtol::Float64,maxit::Int64)

    x   = copy(x_0)
    it  = 0
    res = F.Fnrm(F(x))
    r_0 = res
    println("||F(x)|| = ", res)
    while res > rtol*r_0 + atol && it < maxit
        dx  = -F.jac(x)\F(x)
        x   = x + dx
        it += 1
        res = F.Fnrm(F(x))
        println("||F(x)|| = ", res)
    end
    return x, it
end

#=

Part II: A Parameter ID Problem
         See the notes and video for a description of this application.
=#

################################################################################
# The Objective Function
################################################################################

# The nonlinear least squares objective function
function nonlinear_ls(v_n::Array{Float64,1},
                      t_n::LinRange{Float64},
                        x::Array{Float64,1})
    F = 0.0
    n = length(v_n)
    if n != length(t_n)
        error("v_n and t_n must be the same length")
    end

    it = 1
    while it <= n
        F += (v_n[it] - exp(-x[1]*t_n[it])*(
                             x[3]*sin(x[2]*t_n[it]) +
                             x[4]*cos(x[2]*t_n[it])))^2
        it += 1
    end

    return F
end

################################################################################
# Problem Data
################################################################################

#=
We generate a problem instance
- ex1, ex2, ex3, ex4 represent the "true" unknown parameters
- t is the discretized time interval [0,40] with 800 time steps
- v is are the "true" measurements using the model v(t,x)
=#

# Generate "Data"
ex1 = 0.15
ex2 = 2.0
ex3 = 1.0
ex4 = 3.0

# Generate 800 steps from time 0 to 40
t = LinRange(0, 40, 800)

# v is the vector of time-series data
v = exp.(-ex1.*t).*(ex3*sin.(ex2.*t)+ex4.*cos.(ex2.*t))

theme(:wong)
p1 = plot(t,v,
          titlefont = font(18, "Helvetica"),
          title= "Motion of weight attached to damped spring",
          xlabel="\$t\$",
          ylabel="\$v(t)\$",
          titlefontsize = 24,
          label = false,#["abs_err"],
          bottom_margin=10mm,
          left_margin=14mm,
          right_margin=14mm,
          xtickfontsize=10,
          ytickfontsize=10,
          yguidefontsize=18,
          xguidefontsize=24,
          legend=true,
          xtickfont = font(12, "Helvetica"),
          ytickfont = font(12, "Helvetica"),
          # legendfont = font(14, "Helvetica"),
          lw=3)

savefig(p1,"true_solution_v.pdf")

# To make the problem more interesting, we add noise to the measurments
v = exp.(-ex1.*t).*(ex3*sin.(ex2.*t)+ex4.*cos.(ex2.*t)) + randn(800)/50
plot!(t,v,label = false)

################################################################################
# Defining constructors for the objective
################################################################################

# An objective function type
mutable struct Objective
    data_1::Array{Float64,1}
    data_2::LinRange{Float64}
    obj_fn::Function
end

# Evaluate the objective without  re-passing the fixed data vectors v_n and t_n
function (F::Objective)(x::Array{Float64})
    return F.obj_fn(F.data_1,F.data_2,x)
end

# Create an instance of Objective using the data (v,t)
nls_obj = Objective(v,t,nonlinear_ls)

################################################################################
# Approximating Gradients and Hessians using Finite Differences
################################################################################

# A central difference function for finding partial derivatives
function CD(x::Array{Float64},d::Array{Float64},h::Float64,f::Objective)
    return (f(x+d) - f(x-d))/(2*h)
end

# A function for approximating the gradient using finite differences
function CD_grad(x::Array{Float64},h::Float64,f::Objective)
    n  = length(x)
    it = 1

    gradF = zeros(n)
    while it <= n
        gradF[it] = CD(x,h*(1:n .== it),h,f)
        it += 1
    end

    return gradF
end

# Random point
CD_grad(rand(4),1e-6,nls_obj)
# Value at optimal solution
CD_grad([ex1,ex2,ex3,ex4],1e-6,nls_obj)

# A function for approximating the Hessians using finite differences
function CD_hess(x::Array{Float64,1},h::Float64,f::Objective)
     n = length(x)
     i = 1
     j = 1

    HessF = zeros(n,n)
    while i <= n
        while j <= n
            HessF[i,j] = (CD_grad(x+h*(1:n .== j),h,f)[i] -
                          CD_grad(x-h*(1:n .== j),h,f)[i] +
                          CD_grad(x+h*(1:n .== i),h,f)[j] -
                          CD_grad(x-h*(1:n .== i),h,f)[j])/(4.0*h)
            j+=1
        end
        j  = 1
        i += 1
    end

    return HessF
end

# Random point
CD_hess(rand(4),1e-6,nls_obj)

# Value at optimal solution
CD_hess([ex1,ex2,ex3,ex4],1e-6,nls_obj)

################################################################################
# Defining constructors for gradients and hessians
################################################################################
mutable struct Grad
    f::Objective
    step::Float64
end

function (G::Grad)(x::Array{Float64,1})
    return CD_grad(x,G.step,G.f)
end

mutable struct Hess
    f::Objective
    step::Float64
end

function (H::Hess)(x::Array{Float64,1})
    return CD_hess(x,H.step,H.f)
end

# A composite type for optimization problems
mutable struct OptProb
    obj::Objective
    grd::Grad
    hess::Hess
end

nls_obj  = Objective(v,t,nonlinear_ls)
grad_obj = Grad(nls_obj,1e-8)
hess_obj = Hess(nls_obj,1e-8)

nls_res = OptProb(nls_obj,grad_obj,hess_obj)

# An adaptation of Newton's method for optimization
function newton_nd(F::OptProb,x_0::Array{Float64,1},
                   atol::Float64,rtol::Float64,maxit::Int64)
    x   = copy(x_0)
    it  = 0
    rhs = F.grd(x)
    res = norm(rhs)
    r_0 = res
    println("||F(x)|| = ", res)
    # while res > tol && it < maxit
    while res > rtol*r_0 + atol && it < maxit
        dx  = -F.hess(x)\rhs
        x   = x + dx
        # alpha = 1/(1 + 10*norm(dx)) # Damped version
        # x   = x+alpha*dx
        it += 1
        rhs = F.grd(x)
        res = norm(rhs)
        println("||F(x)|| = ", res)
    end
    println("Number of Iterations = ", it)
    return x, it
end

# Solve the problem!
extn, nit = newton_nd(nls_res,[ex1,ex2,ex3,ex4].+rand(4)/50,1e-8,1e-8,50)

################################################################################
# Plot and compare the approximation
################################################################################
v = exp.(-ex1.*t).*(ex3*sin.(ex2.*t)+ex4.*cos.(ex2.*t))
v_newt = exp.(-extn[1].*t).*(extn[3]*sin.(extn[2].*t)+extn[4].*cos.(extn[2].*t))

p2 = plot(t,v_newt,
          titlefont = font(18, "Helvetica"),
          title= "Reconstructed Motion",
          xlabel="\$t\$",
          ylabel="\$v(t)\$",
          titlefontsize = 24,
          label = false,#["abs_err"],
          bottom_margin=10mm,
          left_margin=14mm,
          right_margin=14mm,
          xtickfontsize=10,
          ytickfontsize=10,
          yguidefontsize=18,
          xguidefontsize=24,
          legend=true,
          xtickfont = font(12, "Helvetica"),
          ytickfont = font(12, "Helvetica"),
          # legendfont = font(14, "Helvetica"),
          lw=3)

savefig(p2,"reconstructed_solution_v.pdf")

err = log.(abs.(v-v_newt))

p3 = plot(t,err,
          titlefont = font(18, "Helvetica"),
          title= "Log error",
          xlabel="\$t\$",
          ylabel="\$\\|v_{true}(t) - v_{newt}(t)|\$",
          titlefontsize = 24,
          label = false,#["abs_err"],
          bottom_margin=10mm,
          left_margin=14mm,
          right_margin=14mm,
          xtickfontsize=10,
          ytickfontsize=10,
          yguidefontsize=18,
          xguidefontsize=24,
          legend=true,
          xtickfont = font(12, "Helvetica"),
          ytickfont = font(12, "Helvetica"),
          # legendfont = font(14, "Helvetica"),
          lw=3)

savefig(p3,"log_error_noisy_data_v.pdf")

p4 = plot(t,(v-v_newt),
          titlefont = font(18, "Helvetica"),
          title= "Absolute error",
          xlabel="\$t\$",
          ylabel="\$\\|v_{true}(t) - v_{newt}(t)|\$",
          titlefontsize = 24,
          label = false,#["abs_err"],
          bottom_margin=10mm,
          left_margin=14mm,
          right_margin=14mm,
          xtickfontsize=10,
          ytickfontsize=10,
          yguidefontsize=18,
          xguidefontsize=24,
          legend=true,
          xtickfont = font(12, "Helvetica"),
          ytickfont = font(12, "Helvetica"),
          # legendfont = font(14, "Helvetica"),
          lw=3)

savefig(p4,"abs_error_noisy_data_v.pdf")

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

#=

Part III: Variational Denoising of Signals
          See the notes and video for a description of this application.
=#

################################################################################
# Generate Residual F(x)
################################################################################
global ts = 1e3      # Penalty of misfit (Try not to use global constants...)

# The nonlinear term in the ODE
function bulk(u::Array{Float64,1},reg::Float64)
    n = length(u)
    F = zeros(n)
    for i in 2:(n-1)
        F[i] = (2*u[i]^3 - 2*u[i]^2 + (1+reg*ts)*u[i])/reg
    end
    return F
end

# This evaluates the full nonlinear residual
function nl_residual(u::Array{Float64,1},
                     g::Array{Float64,1},
                     reg::Float64,
                     L::SparseMatrixCSC)

    return L*u + bulk(u,reg) - [0 ; ts*g ; 0]
end

# Test
# nl_residual(abs.(rand(nnd+2)),abs.(rand(nnd)),reg,L)

################################################################################
# Generate Jacobian
################################################################################
# The nonlinearity is simple enough to directly calcualte and iplement the
# Jacobian
function diff_bulk(u::Array{Float64,1},reg::Float64)
    n = length(u)
    F = zeros(n)
    for i in 2:(n-1)
        F[i] = (6*u[i]^2 - 4*u[i] + (1+ts*reg))/reg
    end
    return F
end

################################################################################
# Defining constructors
################################################################################
# We once again slightly adapt our Newton algorithm by introducing a custom type
mutable struct SigRes
    evalF::Function
    jac::Function
    sig::Array{Float64,1}
    reg::Float64
    L::SparseMatrixCSC
end

# We make Residual a callable type that acts as a function but encapsulates
# data or other important functions.
function (F::SigRes)(x::Array{Float64,1})
    return F.evalF(x,F.sig,F.reg,F.L)
end

# This is the true Hessian for the nonlinear residual
function sig_hess(x::Array{Float64,1},r::Float64,L::SparseMatrixCSC)
    dF  = spdiagm(0 => diff_bulk(x,r))
    return L + dF
end

function newton_dn(F::SigRes,x_0::Array{Float64,1},
                   atol::Float64,rtol::Float64,maxit::Int64)

    x   = copy(x_0)
    it  = 0
    res = norm(F(x))
    r_0 = res
    println("||F(x)|| = ", res)
    while res > rtol*r_0 + atol && it < maxit
        H   = F.jac(x,F.reg,F.L)
        dx  = -H\F(x)
        x   = x + dx
        it += 1
        res = norm(F(x))
        println("||F(x)|| = ", res)
    end
    return x
end

# Solve
nde = 12
nnd = 2^nde          # Number of abscissae
reg = 0.001           # "Smoothing" Parameter for Ginzburg-Landau Function
stp = 1/(2^nde - 1)  # Fixed length of intervals

################################################################################
# Generate a Random Signal Function
################################################################################
function step_function(r_0::Float64,r_1::Float64,s::Float64,x::Float64)
    if x < r_0
        return 0.0
    elseif x >= r_0 && x < r_1
        return s
    else
        return 0.0
    end
end
step_function(0.0,0.25,1.0,-1.0)

function simple_step_func(R::Array{Float64,1},S::Array{Float64,1},x::Float64)
    l = length(R)
    f = 0.0
    for j in 1:(l-1)
        f = f + step_function(R[j],R[j+1],S[j],x)
    end
    return f
end
# Test
R = collect(0:0.1:1)
S = rand((length(R)))
simple_step_func(R,S,0.95)

# Plot
x_grid = LinRange(0,1,nnd)
x_grid_values = collect(x_grid)
signal = zeros(nnd)
clnsig = zeros(nnd)
for i in 1:nnd
    clnsig[i] = simple_step_func(R,S,x_grid_values[i]) # Clean
    signal[i] = clnsig[i] + randn()/100
end

# pp = plot(x_grid,[clnsig,signal],
#           titlefontsize = 11,
#           legend=false,
#           lw=3,
#           )
# display(pp)

# This is the central difference matrix for the term epsilon*u''(x)
L  =  (reg/(stp^2))*spdiagm(-1 =>  -ones(nnd+1),
                             0 => 2*ones(nnd+2),
                             1 =>  -ones(nnd+1))

# This incorporates the homogeneous Newmann boundary conditions into L
L[1,:]     = sparse((1:(nnd+2) .== 1) - (1:(nnd+2) .== 3))
L[nnd+2,:] = sparse((1:(nnd+2) .== (nnd+2)) - (1:(nnd+2) .== nnd))

u_0     = abs.(rand(nnd+2))
sig_obj = SigRes(nl_residual,sig_hess,signal,reg,L)
u_0     = newton_dn(sig_obj,u_0,1e-11,1e-11,10)
plot(u_0[2:(end-1)])

# Compare the results for FD vs. CD
linetypes = [:T :N :DN]
styles = filter((s->begin
                s in Plots.supported_styles()
            end), [:solid, :dash, :dot])
styles = reshape(styles, 1, length(styles))

# The TV regularization is great for denoising but (naturally) affects the
# contrast in the denoised signal. This can be (heuristically!!!) adjusted
# after solving by scaling the solution u, e.g.  3/2 * u gives a decent
# adjustment of contrast.
cntrst = 1.0
# cntrst = 1.5

pp1 = plot(x_grid,[clnsig,signal,cntrst*u_0[2:(nnd+1)]],
          titlefont = font(14, "Helvetica"),
          title= "Noisy vs. Denoised",
          xlabel="\$x\$",
          ylabel="\$u(x)\$",
          titlefontsize = 24,
          label = ["T"  "N"  "DN"],
          line=(1.75, styles),
          bottom_margin=10mm,
          left_margin=14mm,
          right_margin=14mm,
          xtickfontsize=10,
          ytickfontsize=10,
          yguidefontsize=18,
          xguidefontsize=24,
          legend = :bottomleft,
          xtickfont = font(12, "Helvetica"),
          ytickfont = font(12, "Helvetica"),
          legendfont = font(11, "Helvetica")
          )

savefig(pp1,"denoised_signal.pdf")


text:
using LinearAlgebra
using DataFrames
using Test
println("LATHE.PREPROCESS TESTS")
println("......................")
# DECOMPOSITION
using Lathe.preprocess: OneHotEncoder, FloatEncoder, OrdinalEncoder
# ENCODING
df = DataFrame(:A => [5, 10, 15, 20], :B => ["test", "test2", "test2", "test"])
@testset "OneHotEncoder Tests" begin
encoder = OneHotEncoder()
ant = size(df)[2]
data = encoder.predict(df, :A)
len = Set(df[!, :A])
@test size(data)[2] - ant == length(len)

end
@testset "OrdinalEncoder Tests" begin
encoder = OrdinalEncoder(df[!, :B])
encoded = encoder.predict(df[!, :B])
@test encoded[1] == 1
@test encoded[2] == 2
@test encoded[3] == 2
end
@testset "FloatEncoder Tests" begin
  encoder = FloatEncoder()
  encoded = encoder.predict(df[!, :B])
  @test encoded[1] == sum([c = Float64(c) for c in "test"])
end


text:
# Dictionary which returns (and sets) a default value for a requested item not
# already in to the dictionary

# DefaultDictBase is the main class used to in Default*Dicts.
#
# Each related (immutable) Default*Dict class contains a single
# DefaultDictBase object as a member, and delegates almost all
# functions to this object.
#
# The main rationale for doing this instead of using type aliases is
# that this way, we can have actual class names and constructors for
# each of the DefaultDictBase "subclasses", in some sense getting
# around the Julia limitation of not allowing concrete classes to be
# subclassed.
#

struct DefaultDictBase{K,V,F,D} <: AbstractDict{K,V}
    default::F
    d::D
    passkey::Bool

    check_D(D,K,V) = (D <: AbstractDict{K,V}) ||
        throw(ArgumentError("Default dict must be <: AbstractDict{$K,$V}"))

    DefaultDictBase{K,V,F,D}(x::F, kv::AbstractArray{Tuple{K,V}}; passkey=false) where {K,V,F,D} =
        (check_D(D,K,V); new{K,V,F,D}(x, D(kv), passkey))
    DefaultDictBase{K,V,F,D}(x::F, ps::Pair{K,V}...; passkey=false) where {K,V,F,D} =
        (check_D(D,K,V); new{K,V,F,D}(x, D(ps...), passkey))

    DefaultDictBase{K,V,F,D}(x::F, d::D; passkey=d.passkey) where {K,V,F,D<:DefaultDictBase} =
        (check_D(D,K,V); DefaultDictBase(x, d.d; passkey=passkey))
    DefaultDictBase{K,V,F,D}(x::F, d::D = D(); passkey=false) where {K,V,F,D} =
        (check_D(D,K,V); new{K,V,F,D}(x, d, passkey))
end

# Constructors

DefaultDictBase(; kwargs...) = throw(ArgumentError("no default specified"))
DefaultDictBase(k, v; kwargs...) = throw(ArgumentError("no default specified"))

# syntax entry points
DefaultDictBase(default::F; kwargs...) where {F} = DefaultDictBase{Any,Any,F,Dict{Any,Any}}(default; kwargs...)
DefaultDictBase(default::F, kv::AbstractArray{Tuple{K,V}}; kwargs...) where {K,V,F} = DefaultDictBase{K,V,F,Dict{K,V}}(default, kv; kwargs...)
DefaultDictBase(default::F, ps::Pair{K,V}...; kwargs...) where {K,V,F} = DefaultDictBase{K,V,F,Dict{K,V}}(default, ps...; kwargs...)
DefaultDictBase(default::F, d::D; kwargs...) where {F,D<:AbstractDict} = (K=keytype(d); V=valtype(d); DefaultDictBase{K,V,F,D}(default, d; kwargs...))

# Constructor for DefaultDictBase{Int,Float64}(0.0)
DefaultDictBase{K,V}(default::F; kwargs...) where {K,V,F} = DefaultDictBase{K,V,F,Dict{K,V}}(default; kwargs...)

# Functions

# most functions are simply delegated to the wrapped dictionary
@delegate DefaultDictBase.d [ get, haskey, getkey, pop!,
                              iterate, isempty, length ]

# Some functions are delegated, but then need to return the main dictionary
# NOTE: push! is not included below, because the fallback version just
#       calls setindex!
@delegate_return_parent DefaultDictBase.d [ delete!, empty!, setindex!, sizehint! ]

empty(d::DefaultDictBase{K,V,F}) where {K,V,F} = DefaultDictBase{K,V,F}(d.default; passkey=d.passkey)
@deprecate similar(d::DefaultDictBase) empty(d)

getindex(d::DefaultDictBase, key) = get!(d.d, key, d.default)

function getindex(d::DefaultDictBase{K,V,F}, key) where {K,V,F<:Base.Callable}
    if d.passkey
        return get!(d.d, key) do
            d.default(key)
        end
    else
        return get!(d.d, key) do
            d.default()
        end
    end
end

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

# Here begins the actual definition of the DefaultDict and
# DefaultOrderedDict classes.  As noted above, these are simply
# wrappers around a DefaultDictBase object, and delegate all functions
# to that object

for _Dict in [:Dict, :OrderedDict]
    DefaultDict = Symbol("Default"*string(_Dict))
    @eval begin
        struct $DefaultDict{K,V,F} <: AbstractDict{K,V}
            d::DefaultDictBase{K,V,F,$_Dict{K,V}}

            $DefaultDict{K,V,F}(x, ps::Pair{K,V}...; kwargs...) where {K,V,F} =
                new{K,V,F}(DefaultDictBase{K,V,F,$_Dict{K,V}}(x, ps...; kwargs...))
            $DefaultDict{K,V,F}(x, kv::AbstractArray{Tuple{K,V}}; kwargs...) where {K,V,F} =
                new{K,V,F}(DefaultDictBase{K,V,F,$_Dict{K,V}}(x, kv; kwargs...))
            $DefaultDict{K,V,F}(x, d::$DefaultDict) where {K,V,F} = $DefaultDict(x, d.d)
            $DefaultDict{K,V,F}(x, d::$_Dict; kwargs...) where {K,V,F} =
                new{K,V,F}(DefaultDictBase{K,V,F,$_Dict{K,V}}(x, d; kwargs...))
            $DefaultDict{K,V,F}(x; kwargs...) where {K,V,F} =
                new{K,V,F}(DefaultDictBase{K,V,F,$_Dict{K,V}}(x; kwargs...))
        end

        ## Constructors

        $DefaultDict() = throw(ArgumentError("$DefaultDict: no default specified"))
        $DefaultDict(k,v) = throw(ArgumentError("$DefaultDict: no default specified"))

        # syntax entry points
        $DefaultDict(default::F; kwargs...) where {F} = $DefaultDict{Any,Any,F}(default; kwargs...)
        $DefaultDict(default::F, kv::AbstractArray{Tuple{K,V}}; kwargs...) where {K,V,F} = $DefaultDict{K,V,F}(default, kv; kwargs...)
        $DefaultDict(default::F, ps::Pair{K,V}...; kwargs...) where {K,V,F} = $DefaultDict{K,V,F}(default, ps...; kwargs...)

        $DefaultDict(default::F, d::AbstractDict; kwargs...) where {F} = ((K,V)= (Base.keytype(d), Base.valtype(d)); $DefaultDict{K,V,F}(default, $_Dict(d); kwargs...))

        # Constructor syntax: DefaultDictBase{Int,Float64}(default)
        $DefaultDict{K,V}(; kwargs...) where {K,V} = throw(ArgumentError("$DefaultDict: no default specified"))
        $DefaultDict{K,V}(default::F; kwargs...) where {K,V,F} = $DefaultDict{K,V,F}(default; kwargs...)

        ## Functions

        # Most functions are simply delegated to the wrapped DefaultDictBase object
        @delegate $DefaultDict.d [ getindex, get, get!, haskey,
                                   getkey, pop!, iterate,
                                   isempty, length ]

        # Some functions are delegated, but then need to return the main dictionary
        # NOTE: push! is not included below, because the fallback version just
        #       calls setindex!
        @delegate_return_parent $DefaultDict.d [ delete!, empty!, setindex!, sizehint! ]

        # NOTE: The second and third definition of push! below are only
        # necessary for disambiguating with the fourth, fifth, and sixth
        # definitions of push! below.
        # If these are removed, the second and third definitions can be
        # removed as well.
        push!(d::$DefaultDict, p::Pair) = (setindex!(d.d, p.second, p.first); d)
        push!(d::$DefaultDict, p::Pair, q::Pair) = push!(push!(d, p), q)
        push!(d::$DefaultDict, p::Pair, q::Pair, r::Pair...) = push!(push!(push!(d, p), q), r...)

        push!(d::$DefaultDict, p) = (setindex!(d.d, p[2], p[1]); d)
        push!(d::$DefaultDict, p, q) = push!(push!(d, p), q)
        push!(d::$DefaultDict, p, q, r...) = push!(push!(push!(d, p), q), r...)

        empty(d::$DefaultDict{K,V,F}) where {K,V,F} = $DefaultDict{K,V,F}(d.d.default)
        in(key, v::Base.KeySet{K,T}) where {K,T<:$DefaultDict{K}} = key in keys(v.dict.d.d)

        @deprecate similar(d::$DefaultDict) empty(d)
    end
end

isordered(::Type{T}) where {T<:DefaultOrderedDict} = true

## This should be uncommented to provide a DefaultSortedDict

# struct DefaultSortedDict{K,V,F} <: AbstractDict{K,V}
#     d::DefaultDictBase{K,V,F,SortedDict{K,V}}

#     DefaultSortedDict(x, kv::AbstractArray{(K,V)}) = new(DefaultDictBase{K,V,F,SortedDict{K,V}}(x, kv))
#     DefaultSortedDict(x, d::DefaultSortedDict) = DefaultSortedDict(x, d.d)
#     DefaultSortedDict(x, d::SortedDict) = new(DefaultDictBase{K,V,F,SortedDict{K,V}}(x, d))
#     DefaultSortedDict(x) = new(DefaultDictBase{K,V,F,SortedDict{K,V}}(x))
#     DefaultSortedDict(x, ks, vs) = new(DefaultDictBase{K,V,F,SortedDict{K,V}}(x,ks,vs))
# end

## Constructors

# DefaultSortedDict() = throw(ArgumentError("DefaultSortedDict: no default specified"))
# DefaultSortedDict(k,v) = throw(ArgumentError("DefaultSortedDict: no default specified"))

# # TODO: these mimic similar Dict constructors, but may not be needed
# DefaultSortedDict{K,V,F}(default::F, ks::AbstractArray{K}, vs::AbstractArray{V}) = DefaultSortedDict{K,V,F}(default,ks,vs)
# DefaultSortedDict{F}(default::F,ks,vs) = DefaultSortedDict{Any,Any,F}(default, ks, vs)

# # syntax entry points
# DefaultSortedDict{F}(default::F) = DefaultSortedDict{Any,Any,F}(default)
# DefaultSortedDict{K,V,F}(::Type{K}, ::Type{V}, default::F) = DefaultSortedDict{K,V,F}(default)
# DefaultSortedDict{K,V,F}(default::F, kv::AbstractArray{(K,V)}) = DefaultSortedDict{K,V,F}(default, kv)

# DefaultSortedDict{F}(default::F, d::AbstractDict) = ((K,V)=eltype(d); DefaultSortedDict{K,V,F}(default, SortedDict(d)))

## Functions

## Most functions are simply delegated to the wrapped DefaultDictBase object

# @delegate DefaultSortedDict.d [ getindex, get, get!, haskey,
#                                 getkey, pop!, start, next,
#                                 done, isempty, length]
# @delegate_return_parent DefaultSortedDict.d [ delete!, empty!, setindex!, sizehint! ]

# similar{K,V,F}(d::DefaultSortedDict{K,V,F}) = DefaultSortedDict{K,V,F}(d.d.default)
# in{T<:DefaultSortedDict}(key, v::Base.KeySet{T}) = key in keys(v.dict.d.d)


text:
using GeometricIntegrators

set_config(:tab_compensated_summation, true)

const out_dir = "../../docs/src/integrators/projection_methods/guiding_center_4d/"
const run_id  = "fast_barely_trapped_vprk_comp"

include("guiding_center_4d_fast_barely_trapped.jl")
include("../guiding_center_4d_vprk.jl")


text:
# Test rosenbrock pdf gradients
# Authors: Ali Siahkoohi, alisk@gatech.edu
# Date: September 2020

using DrWatson
@quickactivate :FastApproximateInference
import Pkg; Pkg.instantiate()

using LinearAlgebra
using Test

# Input
n_samples = 1

# Rosenbrock distribution
RB_dist = RosenbrockDistribution(0f0, 5f-1, ones(Float32, 2, 1))

# samples
X0 = rand(RB_dist, n_samples)
X = rand(RB_dist, n_samples)
dX = X - X0

function rosenbrock_nlpdf(RB_dist, x)
    return nl_pdf(x, RB_dist)[1], nl_pdf_grad(x, RB_dist)
end

# Gradient test
f0, ΔX = rosenbrock_nlpdf(RB_dist, X0)
h = 0.1f0
maxiter = 5
err1 = zeros(Float32, maxiter)
err2 = zeros(Float32, maxiter)

print("\nGradient test\n")
for j=1:maxiter
    f = rosenbrock_nlpdf(RB_dist, X0 + h*dX)[1]
    err1[j] = abs(f - f0)
    err2[j] = abs(f - f0 - h*dot(dX, ΔX))
    print(err1[j], "; ", err2[j], "\n")
    global h = h/2f0
end

@test isapprox(err1[end] / (err1[1]/2^(maxiter-1)), 1f0; atol=1f1)
@test isapprox(err2[end] / (err2[1]/4^(maxiter-1)), 1f0; atol=1f1)


text:
module PointerArithmetic

    @warn "Pointer and ArraySlice are unsafe! \n" *
          "In order to allow you working with pointer arithmetics, \n" *
          "they don't have bounds checking, and it can lead to OS \n" *
          "or application malfunctioning if used wrongly. \n" *
          "Make sure you know what you're doing before using them." maxlog = 1

    using Reexport

    include("Misc.jl")
    using .Misc

    include("Pointer.jl")
    @reexport using .Pointers

    include("ArraySlice.jl")
    @reexport using .ArraySlices

    include("Slice.jl")
    @reexport using .Slices

end


text:
#=
This file is auto-generated. Do not edit.
=#
"""
    mutable struct RenewableFix <: RenewableGen
        name::String
        available::Bool
        bus::Bus
        activepower::Float64
        reactivepower::Float64
        tech::TechRenewable
        services::Vector{Service}
        ext::Dict{String, Any}
        forecasts::InfrastructureSystems.Forecasts
        internal::InfrastructureSystemsInternal
    end

Data Structure for fixed renewable generation technologies.

# Arguments
- `name::String`
- `available::Bool`
- `bus::Bus`
- `activepower::Float64`
- `reactivepower::Float64`
- `tech::TechRenewable`
- `services::Vector{Service}`: Services that this device contributes to
- `ext::Dict{String, Any}`
- `forecasts::InfrastructureSystems.Forecasts`: internal forecast storage
- `internal::InfrastructureSystemsInternal`: power system internal reference, do not modify
"""
mutable struct RenewableFix <: RenewableGen
    name::String
    available::Bool
    bus::Bus
    activepower::Float64
    reactivepower::Float64
    tech::TechRenewable
    "Services that this device contributes to"
    services::Vector{Service}
    ext::Dict{String, Any}
    "internal forecast storage"
    forecasts::InfrastructureSystems.Forecasts
    "power system internal reference, do not modify"
    internal::InfrastructureSystemsInternal
end

function RenewableFix(name, available, bus, activepower, reactivepower, tech, services=Device[], ext=Dict{String, Any}(), forecasts=InfrastructureSystems.Forecasts(), )
    RenewableFix(name, available, bus, activepower, reactivepower, tech, services, ext, forecasts, InfrastructureSystemsInternal(), )
end

function RenewableFix(; name, available, bus, activepower, reactivepower, tech, services=Device[], ext=Dict{String, Any}(), forecasts=InfrastructureSystems.Forecasts(), )
    RenewableFix(name, available, bus, activepower, reactivepower, tech, services, ext, forecasts, )
end

# Constructor for demo purposes; non-functional.
function RenewableFix(::Nothing)
    RenewableFix(;
        name="init",
        available=false,
        bus=Bus(nothing),
        activepower=0.0,
        reactivepower=0.0,
        tech=TechRenewable(nothing),
        services=Device[],
        ext=Dict{String, Any}(),
        forecasts=InfrastructureSystems.Forecasts(),
    )
end

InfrastructureSystems.get_name(value::RenewableFix) = value.name
"""Get RenewableFix available."""
get_available(value::RenewableFix) = value.available
"""Get RenewableFix bus."""
get_bus(value::RenewableFix) = value.bus
"""Get RenewableFix activepower."""
get_activepower(value::RenewableFix) = value.activepower
"""Get RenewableFix reactivepower."""
get_reactivepower(value::RenewableFix) = value.reactivepower
"""Get RenewableFix tech."""
get_tech(value::RenewableFix) = value.tech
"""Get RenewableFix services."""
get_services(value::RenewableFix) = value.services
"""Get RenewableFix ext."""
get_ext(value::RenewableFix) = value.ext

InfrastructureSystems.get_forecasts(value::RenewableFix) = value.forecasts
"""Get RenewableFix internal."""
get_internal(value::RenewableFix) = value.internal