text:
using FantasyFootball
using Gadfly, Color
using DataFrames, DataFramesMeta
using StatsBase

# include stat functions that can work around NAs
include("stat.jl")

gradient = Scale.ContinuousColorScale(Scale.lab_gradient(color("white"), RGB(0.5, 0, 0), color("red")))

# target dir should be checked into git, but empty.
dir = joinpath(dirname(@__FILE__), "..", "graphs")
# TODO - clean target directory here?

# General experiments

# General setup
s10 = score(Season(2010))
s11 = score(Season(2011))
s12 = score(Season(2012))
s13 = score(Season(2013))
s14 = score(Season(2014))

seasons = [s10, s11, s12, s13, s14]

off = [[@transform(getfield(seasons[i].off, f), pos=string(f), season=string(2009 + i)) for f in names(seasons[i].off)] for i in 1:length(seasons)]

fields = [:team_abbr, :score]
plucked = [
  [@transform(getfield(s14.off, f)[fields], pos = string(f)) for f in names(s14.off)],
  @transform(s14.def.def[fields], pos = "def/st"),
  @transform(s14.k.k[fields], pos = "k")
]

df = vcat(plucked...)
#
# pos_tables = [nu
# for season in seasons:
#   off = @

# Point density by position
# pd_pos = plot(df, x = :score, color = :pos, Geom.density, Guide.title("Point production density by position, 2014"))
# draw(SVG(joinpath(dir, "pos_point_density.svg"), 30cm, 12cm), pd_pos)
# hist_pos = plot(df, x = :score, color = :pos, Geom.histogram, Guide.title("Point production by position"))
# draw(SVG(joinpath(dir, "pos_point_hist.svg"), 30cm, 12cm), hist_pos)

# # "Rounds" pos density
# rounds = @transform(sort(df, cols = [:score], rev = true), round = [div(1, x) + 1 for x in 1:size(df)[1]])
# pd_round = plot(@where(rounds, :round .< 10), x = :round, y = :pos, Geom.histogram2d, Scale.x_discrete, gradient, Guide.title("Position distribution of top hundred picks by draft round"))
# #draw(SVG(joinpath(dir, "top_positions_by_rounds.svg"), 12cm, 6cm), pd_round)

# Names for the point columns
pt_cols = [symbol("pts$(s)") for s in 2010:2014]

# Helpers to extract score info by position and join the year tables
# e.g.
# qbs = [sort(s.off.qb), cols = [:score], rev = true)[[:Name, :score]] for s in seasons]
# all_qb = reduce((a, b) -> join(a, b, on = :Name, kind = :outer), qbs)
join_pos(pos, j) = reduce((a, b) -> join(a, b, on = j, kind = :outer), pos)
get_tables(pos, j) = [sort(getfield(s.off, pos), cols = [:score], rev = true)[[j, :score]] for s in seasons]
get_tables(p1, p2, j) = [sort(getfield(getfield(s, p1), p2), cols = [:score], rev = true)[[j, :score]] for s in seasons]
function get_all_years(p1, p2, j=:Name)
    table = join_pos(get_tables(p1, p2), j)
    names!(table, [j, pt_cols])
    @transform(table, std2010=nazscore(:pts2010), std2011=nazscore(:pts2011), std2012=nazscore(:pts2012), std2013=nazscore(:pts2013), std2014=nazscore(:pts2014))
end

wrs = get_all_years(:off, :wr)
rbs = get_all_years(:off, :rb)
tes = get_all_years(:off, :te)
qbs = get_all_years(:off, :qb)
def = get_all_years(:def, :def, :team_abbr)

std_cols = [:std2010, :std2011, :std2012, :std2013, :std2014]
function best_of_2014(table, n=100)
    top = sort(@where(table, !isna(:std2014)), cols=[sym], rev=true)[[:Name, std_cols]]
    top_lng = @where(melt(top[1:n,:], :Name), !isna(:value))
end

import Base.max
function max{T<:Real}(cols::DataArray{T}...)
    Z = DataArray(T, size(cols[1]))
    for i in 1:length(cols[1])
        @inbounds Z[i] = max([isna(col[i]) ? -Inf : col[i] for col in cols]...)
    end
    Z
end

function best_ever(table, n)
    tbl = @transform(table[[:Name, std_cols]], max = max(:std2010, :std2011, :std2012, :std2013, :std2014))
    best = sort(tbl, cols=[:max], rev=true)[1:n,:]
    top_melt = melt(best[[:Name, std_cols]], :Name)
    sort(@where(top_melt, !isna(:value)), cols=[:variable, :value])
end

_top_line(df) = plot(df, x=:variable, y=:value, color=:Name, Geom.line)
_top_rect(df) = plot(df, x=:variable, y=:Name, color=:value, Geom.rectbin)

top2014_line(table, n=10) = _top_line(best_of_2014(table, n))
top2014_rect(table, n=20) = _top_rect(best_of_2014(table, n))
best_line(table, n=20) = _top_line(best_ever(table, n))
best_rect(table, n=20) = _top_rect(best_ever(table, n))

max_stars_repo_name: jnu/FantasyFootball

max_issues_repo_name: jnu/FantasyFootball

max_forks_repo_name: jnu/FantasyFootball


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

using Markdown
using InteractiveUtils

# ╔═╡ c7c324de-a354-11eb-287e-2d1d6a6ab5e0
begin
    import Pkg

    Pkg.activate(mktempdir())
    Pkg.add([
        Pkg.PackageSpec(name="Plots", version="1"),
        Pkg.PackageSpec(name="PlutoUI", version="0.7"),
        Pkg.PackageSpec(name="StatsBase", version="0.33"),
    ])

    using Plots, PlutoUI, StatsBase, Statistics
end

# ╔═╡ ef0d7c89-6017-4f42-8a85-db24dafb67e1
PlutoUI.TableOfContents(aside=true)

# ╔═╡ fb7c54c6-79d8-41f9-bb51-1bbb3193588c
md"""
# Modelling a simple Pendulum
"""

# ╔═╡ 4df74a52-c618-4e04-b44b-0d610c3b58a0
md"""
Modelling the motion of a simple pendulum consisting of a stiff, inextensible rod of length l with a bob of mass m. It's position can be determined by the angle $\theta$ that it makes with the central reference line.

Let the gravtitational acceleration be denoted by $g$.

Restorative force on the pendulum is given by,

$$F_R = - m g sin(\theta)$$

This is equal to mass times acceleration a,

$$F = m a$$

$$m a = -m g sin(\theta)$$

$$a = - g sin(\theta)$$

Acceleration is second derivative of displacement s,

$$a = \frac{d^2 s}{dt^2}$$

where

$$s = l \theta$$

Hence

$$a = l \frac{d^2 \theta}{dt^2}$$

Let the velocity be denoted by v then it can be represented as,

$$l \frac{d^2 \theta}{dt^2} = - g sin(\theta)$$

$$\frac{d^2 \theta}{dt^2} = - \frac{g}{l} sin(\theta)$$

or

$$\ddot{\theta} = - \frac{g}{l} sin(\theta)$$
"""

# ╔═╡ d1af15e4-2beb-4a7f-9388-d78523c433a7
md"""
## Euler Method to find pendulum position
"""

# ╔═╡ 3ea4dc17-4868-404f-a575-b97a79badbb1
md"""
Let the initial position of the pendulum be $\theta_0$ and we want to find to the position of the pendulum in terms of the angle formed with the vertical reference line on $k^{th}$ day as $\theta_k$.

The position on the pendulum on day k+1, can be found using the position

Using the Euler method to discreetize the equation above,

$$\ddot{\theta} = - \frac{g}{l} sin(\theta)$$

Approximate the derivative using a small (but not *too* small) time step $\delta t$:

$$\ddot{\theta} \simeq \frac{\dot{\theta}(t + \delta t) - \dot{\theta}(t)}{\delta t}$$

giving

$$\dot{\theta}(t + \delta t) \simeq \dot{\theta}(t) + \delta t \, \ddot{\theta}$$

and

$$\theta(t + \delta t) \simeq \theta(t) + \delta t \, \dot{\theta}(t)$$
"""

# ╔═╡ 0e051bc8-e67e-4ca5-8bb4-1f238bbe3f7b
md"""
Initial value conditions for the pendulum model are,

Initial position at t=0,

$$\theta(t=0) = \theta_0$$

$$\dot{\theta}(t=0) = 0$$
"""

# ╔═╡ 125194cb-f28b-4829-ac5b-484399b3a04c
begin
	θ = Vector{Float64}()
	me_t = Vector{Float64}()
end

# ╔═╡ 203db499-963a-4fa3-957f-b905424d96d9
let
	m = 1
	g = 9.8
	l = 10

	T = 20
	δt = 0.1
	t = 0

	θ_0 = π/12
	θ̇_0 = 0
	θ̈_0 = - (g / l) * sin(θ_0)

	time_range = 0+δt:δt:T

	axis_lim = 10

	append!(θ, θ_0)

	x = l * sin(θ_0)
	y = - l * cos(θ_0)

	PE = m * g * l * cos(θ_0)
	KE = 0.5 * m * l * l * (θ̇_0 * θ̇_0)

	ME = PE + KE

	append!(me_t, ME)

	plot([0,x], [0,y],size=(400,300),xlim=(-2*axis_lim,2*axis_lim),ylim=(-1.5*axis_lim,axis_lim),markersize = 8, markershape = :circle,label ="",axis = []);

	@gif for t in time_range
		str = string("Time = ", round(t, digits = 2), " sec");

		θ_t = θ_0 + (δt * θ̇_0)
		append!(θ, θ_t)

		x = l * sin(θ_t)
		y = - l * cos(θ_t)

		plot!([0,x], [0,y],size=(400,300),xlim=(-2*axis_lim,2*axis_lim),ylim=(-1.5*axis_lim,axis_lim),markersize = 8, markershape = :circle,label ="",axis = [], title = str, title_location = :left);

		θ̇_0 = θ̇_0 + (δt * θ̈_0)
		θ_0 = θ_t
		θ̈_0 = - (g / l) * sin(θ_0)

		PE = m * g * l * cos(θ_0)
		KE = 0.5 * m * l * l * (θ̇_0 * θ̇_0)

		ME = PE + KE
		append!(me_t, ME)

	end
end

# ╔═╡ cb6cebbb-803b-4466-8ae9-5368773e669a
plot(θ, alpha=0.5, ms=2, label="Angle of pendulum with the vertical reference line")

# ╔═╡ 7677fbec-78a6-4486-b10f-4466dffc2c83
plot(me_t, alpha=0.5, ms=2, label="Energy in the simple pendulum")

# ╔═╡ 445440ba-262f-48e2-8c70-5f765beda956
θ

# ╔═╡ 1ac63819-9e52-4e77-b325-109775d4f04d
let
	T = 17632
	δt = 4
	time_range = 0+δt:δt:T
	x = Vector{Float64}()
	append!(x, 0)
	for t in time_range
		append!(x, t)
	end
	[size(time_range, 1)+1 size(x, 1)]
	0+δt:δt:T
end

# ╔═╡ 3bb8acbf-a700-414f-ac9b-f53847b6fc77
begin
	T1 = 100
	δt1 = 0.1
	time_range1 = 0+δt1:δt1:T1
	x1 = Vector{Float64}()
	for t in time_range1
		append!(x1,t)
		if t>3
			append!(x1, t)
		end
	end
	x1
end

# ╔═╡ Cell order:
# ╠═c7c324de-a354-11eb-287e-2d1d6a6ab5e0
# ╠═ef0d7c89-6017-4f42-8a85-db24dafb67e1
# ╟─fb7c54c6-79d8-41f9-bb51-1bbb3193588c
# ╟─4df74a52-c618-4e04-b44b-0d610c3b58a0
# ╟─d1af15e4-2beb-4a7f-9388-d78523c433a7
# ╟─3ea4dc17-4868-404f-a575-b97a79badbb1
# ╟─0e051bc8-e67e-4ca5-8bb4-1f238bbe3f7b
# ╟─125194cb-f28b-4829-ac5b-484399b3a04c
# ╟─203db499-963a-4fa3-957f-b905424d96d9
# ╠═cb6cebbb-803b-4466-8ae9-5368773e669a
# ╠═7677fbec-78a6-4486-b10f-4466dffc2c83
# ╠═445440ba-262f-48e2-8c70-5f765beda956
# ╟─1ac63819-9e52-4e77-b325-109775d4f04d
# ╟─3bb8acbf-a700-414f-ac9b-f53847b6fc77

ext: jl


text:
include("wabt.jl")
include("binaryen.jl")


text:
struct JSONLD_DOI <: JSONLD end

function mainpage_url(repo::JSONLD_DOI, dataname)
    if match_doi(dataname) != nothing
        url = joinpath("https://data.datacite.org/", match_doi(dataname))
        resp = HTTP.get(url, ["Accept"=>"application/vnd.schemaorg.ld+json"]; forwardheaders=true)
        json = resp.body |> String |> JSON.parse
        return json, dataname
    else
        throw(GeneratorError(repo))
    end
end


text:
# Julia wrapper for header: /Users/gnimuc/.julia/dev/LibCSFML/gen/../deps/usr/include/SFML/Window.h
# Automatically generated using Clang.jl

# Julia wrapper for header: /Users/gnimuc/.julia/dev/LibCSFML/gen/../deps/usr/include/SFML/Window/Clipboard.h
# Automatically generated using Clang.jl

function sfClipboard_getString()
    ccall((:sfClipboard_getString, libcsfml_window), Cstring, ())
end

function sfClipboard_getUnicodeString()
    ccall((:sfClipboard_getUnicodeString, libcsfml_window), Ptr{sfUint32}, ())
end

function sfClipboard_setString(text)
    ccall((:sfClipboard_setString, libcsfml_window), Cvoid, (Cstring,), text)
end

function sfClipboard_setUnicodeString(text)
    ccall((:sfClipboard_setUnicodeString, libcsfml_window), Cvoid, (Ptr{sfUint32},), text)
end
# Julia wrapper for header: /Users/gnimuc/.julia/dev/LibCSFML/gen/../deps/usr/include/SFML/Window/Context.h
# Automatically generated using Clang.jl

function sfContext_create()
    ccall((:sfContext_create, libcsfml_window), Ptr{sfContext}, ())
end

function sfContext_destroy(context)
    ccall((:sfContext_destroy, libcsfml_window), Cvoid, (Ptr{sfContext},), context)
end

function sfContext_setActive(context, active)
    ccall((:sfContext_setActive, libcsfml_window), sfBool, (Ptr{sfContext}, sfBool), context, active)
end

function sfContext_getSettings(context)
    ccall((:sfContext_getSettings, libcsfml_window), sfContextSettings, (Ptr{sfContext},), context)
end

function sfContext_getActiveContextId()
    ccall((:sfContext_getActiveContextId, libcsfml_window), sfUint64, ())
end
# Julia wrapper for header: /Users/gnimuc/.julia/dev/LibCSFML/gen/../deps/usr/include/SFML/Window/Cursor.h
# Automatically generated using Clang.jl

function sfCursor_createFromPixels(pixels, size, hotspot)
    ccall((:sfCursor_createFromPixels, libcsfml_window), Ptr{sfCursor}, (Ptr{sfUint8}, sfVector2u, sfVector2u), pixels, size, hotspot)
end

function sfCursor_createFromSystem(type)
    ccall((:sfCursor_createFromSystem, libcsfml_window), Ptr{sfCursor}, (sfCursorType,), type)
end

function sfCursor_destroy(cursor)
    ccall((:sfCursor_destroy, libcsfml_window), Cvoid, (Ptr{sfCursor},), cursor)
end
# Julia wrapper for header: /Users/gnimuc/.julia/dev/LibCSFML/gen/../deps/usr/include/SFML/Window/Event.h
# Automatically generated using Clang.jl

# Julia wrapper for header: /Users/gnimuc/.julia/dev/LibCSFML/gen/../deps/usr/include/SFML/Window/Export.h
# Automatically generated using Clang.jl

# Julia wrapper for header: /Users/gnimuc/.julia/dev/LibCSFML/gen/../deps/usr/include/SFML/Window/Joystick.h
# Automatically generated using Clang.jl

function sfJoystick_isConnected(joystick)
    ccall((:sfJoystick_isConnected, libcsfml_window), sfBool, (UInt32,), joystick)
end

function sfJoystick_getButtonCount(joystick)
    ccall((:sfJoystick_getButtonCount, libcsfml_window), UInt32, (UInt32,), joystick)
end

function sfJoystick_hasAxis(joystick, axis)
    ccall((:sfJoystick_hasAxis, libcsfml_window), sfBool, (UInt32, sfJoystickAxis), joystick, axis)
end

function sfJoystick_isButtonPressed(joystick, button)
    ccall((:sfJoystick_isButtonPressed, libcsfml_window), sfBool, (UInt32, UInt32), joystick, button)
end

function sfJoystick_getAxisPosition(joystick, axis)
    ccall((:sfJoystick_getAxisPosition, libcsfml_window), Cfloat, (UInt32, sfJoystickAxis), joystick, axis)
end

function sfJoystick_getIdentification(joystick)
    ccall((:sfJoystick_getIdentification, libcsfml_window), sfJoystickIdentification, (UInt32,), joystick)
end

function sfJoystick_update()
    ccall((:sfJoystick_update, libcsfml_window), Cvoid, ())
end
# Julia wrapper for header: /Users/gnimuc/.julia/dev/LibCSFML/gen/../deps/usr/include/SFML/Window/JoystickIdentification.h
# Automatically generated using Clang.jl

# Julia wrapper for header: /Users/gnimuc/.julia/dev/LibCSFML/gen/../deps/usr/include/SFML/Window/Keyboard.h
# Automatically generated using Clang.jl

function sfKeyboard_isKeyPressed(key)
    ccall((:sfKeyboard_isKeyPressed, libcsfml_window), sfBool, (sfKeyCode,), key)
end

function sfKeyboard_setVirtualKeyboardVisible(visible)
    ccall((:sfKeyboard_setVirtualKeyboardVisible, libcsfml_window), Cvoid, (sfBool,), visible)
end
# Julia wrapper for header: /Users/gnimuc/.julia/dev/LibCSFML/gen/../deps/usr/include/SFML/Window/Mouse.h
# Automatically generated using Clang.jl

function sfMouse_isButtonPressed(button)
    ccall((:sfMouse_isButtonPressed, libcsfml_window), sfBool, (sfMouseButton,), button)
end

function sfMouse_getPosition(relativeTo)
    ccall((:sfMouse_getPosition, libcsfml_window), sfVector2i, (Ptr{sfWindow},), relativeTo)
end

function sfMouse_setPosition(position, relativeTo)
    ccall((:sfMouse_setPosition, libcsfml_window), Cvoid, (sfVector2i, Ptr{sfWindow}), position, relativeTo)
end
# Julia wrapper for header: /Users/gnimuc/.julia/dev/LibCSFML/gen/../deps/usr/include/SFML/Window/Sensor.h
# Automatically generated using Clang.jl

function sfSensor_isAvailable(sensor)
    ccall((:sfSensor_isAvailable, libcsfml_window), sfBool, (sfSensorType,), sensor)
end

function sfSensor_setEnabled(sensor, enabled)
    ccall((:sfSensor_setEnabled, libcsfml_window), Cvoid, (sfSensorType, sfBool), sensor, enabled)
end

function sfSensor_getValue(sensor)
    ccall((:sfSensor_getValue, libcsfml_window), sfVector3f, (sfSensorType,), sensor)
end
# Julia wrapper for header: /Users/gnimuc/.julia/dev/LibCSFML/gen/../deps/usr/include/SFML/Window/Touch.h
# Automatically generated using Clang.jl

function sfTouch_isDown(finger)
    ccall((:sfTouch_isDown, libcsfml_window), sfBool, (UInt32,), finger)
end

function sfTouch_getPosition(finger, relativeTo)
    ccall((:sfTouch_getPosition, libcsfml_window), sfVector2i, (UInt32, Ptr{sfWindow}), finger, relativeTo)
end
# Julia wrapper for header: /Users/gnimuc/.julia/dev/LibCSFML/gen/../deps/usr/include/SFML/Window/Types.h
# Automatically generated using Clang.jl

# Julia wrapper for header: /Users/gnimuc/.julia/dev/LibCSFML/gen/../deps/usr/include/SFML/Window/VideoMode.h
# Automatically generated using Clang.jl

function sfVideoMode_getDesktopMode()
    ccall((:sfVideoMode_getDesktopMode, libcsfml_window), sfVideoMode, ())
end

function sfVideoMode_getFullscreenModes(count)
    ccall((:sfVideoMode_getFullscreenModes, libcsfml_window), Ptr{sfVideoMode}, (Ptr{Csize_t},), count)
end

function sfVideoMode_isValid(mode)
    ccall((:sfVideoMode_isValid, libcsfml_window), sfBool, (sfVideoMode,), mode)
end
# Julia wrapper for header: /Users/gnimuc/.julia/dev/LibCSFML/gen/../deps/usr/include/SFML/Window/Window.h
# Automatically generated using Clang.jl

function sfWindow_create(mode, title, style, settings)
    ccall((:sfWindow_create, libcsfml_window), Ptr{sfWindow}, (sfVideoMode, Cstring, sfUint32, Ptr{sfContextSettings}), mode, title, style, settings)
end

function sfWindow_createUnicode(mode, title, style, settings)
    ccall((:sfWindow_createUnicode, libcsfml_window), Ptr{sfWindow}, (sfVideoMode, Ptr{sfUint32}, sfUint32, Ptr{sfContextSettings}), mode, title, style, settings)
end

function sfWindow_createFromHandle(handle, settings)
    ccall((:sfWindow_createFromHandle, libcsfml_window), Ptr{sfWindow}, (sfWindowHandle, Ptr{sfContextSettings}), handle, settings)
end

function sfWindow_destroy(window)
    ccall((:sfWindow_destroy, libcsfml_window), Cvoid, (Ptr{sfWindow},), window)
end

function sfWindow_close(window)
    ccall((:sfWindow_close, libcsfml_window), Cvoid, (Ptr{sfWindow},), window)
end

function sfWindow_isOpen(window)
    ccall((:sfWindow_isOpen, libcsfml_window), sfBool, (Ptr{sfWindow},), window)
end

function sfWindow_getSettings(window)
    ccall((:sfWindow_getSettings, libcsfml_window), sfContextSettings, (Ptr{sfWindow},), window)
end

function sfWindow_pollEvent(window, event)
    ccall((:sfWindow_pollEvent, libcsfml_window), sfBool, (Ptr{sfWindow}, Ptr{sfEvent}), window, event)
end

function sfWindow_waitEvent(window, event)
    ccall((:sfWindow_waitEvent, libcsfml_window), sfBool, (Ptr{sfWindow}, Ptr{sfEvent}), window, event)
end

function sfWindow_getPosition(window)
    ccall((:sfWindow_getPosition, libcsfml_window), sfVector2i, (Ptr{sfWindow},), window)
end

function sfWindow_setPosition(window, position)
    ccall((:sfWindow_setPosition, libcsfml_window), Cvoid, (Ptr{sfWindow}, sfVector2i), window, position)
end

function sfWindow_getSize(window)
    ccall((:sfWindow_getSize, libcsfml_window), sfVector2u, (Ptr{sfWindow},), window)
end

function sfWindow_setSize(window, size)
    ccall((:sfWindow_setSize, libcsfml_window), Cvoid, (Ptr{sfWindow}, sfVector2u), window, size)
end

function sfWindow_setTitle(window, title)
    ccall((:sfWindow_setTitle, libcsfml_window), Cvoid, (Ptr{sfWindow}, Cstring), window, title)
end

function sfWindow_setUnicodeTitle(window, title)
    ccall((:sfWindow_setUnicodeTitle, libcsfml_window), Cvoid, (Ptr{sfWindow}, Ptr{sfUint32}), window, title)
end

function sfWindow_setIcon(window, width, height, pixels)
    ccall((:sfWindow_setIcon, libcsfml_window), Cvoid, (Ptr{sfWindow}, UInt32, UInt32, Ptr{sfUint8}), window, width, height, pixels)
end

function sfWindow_setVisible(window, visible)
    ccall((:sfWindow_setVisible, libcsfml_window), Cvoid, (Ptr{sfWindow}, sfBool), window, visible)
end

function sfWindow_setVerticalSyncEnabled(window, enabled)
    ccall((:sfWindow_setVerticalSyncEnabled, libcsfml_window), Cvoid, (Ptr{sfWindow}, sfBool), window, enabled)
end

function sfWindow_setMouseCursorVisible(window, visible)
    ccall((:sfWindow_setMouseCursorVisible, libcsfml_window), Cvoid, (Ptr{sfWindow}, sfBool), window, visible)
end

function sfWindow_setMouseCursorGrabbed(window, grabbed)
    ccall((:sfWindow_setMouseCursorGrabbed, libcsfml_window), Cvoid, (Ptr{sfWindow}, sfBool), window, grabbed)
end

function sfWindow_setMouseCursor(window, cursor)
    ccall((:sfWindow_setMouseCursor, libcsfml_window), Cvoid, (Ptr{sfWindow}, Ptr{sfCursor}), window, cursor)
end

function sfWindow_setKeyRepeatEnabled(window, enabled)
    ccall((:sfWindow_setKeyRepeatEnabled, libcsfml_window), Cvoid, (Ptr{sfWindow}, sfBool), window, enabled)
end

function sfWindow_setFramerateLimit(window, limit)
    ccall((:sfWindow_setFramerateLimit, libcsfml_window), Cvoid, (Ptr{sfWindow}, UInt32), window, limit)
end

function sfWindow_setJoystickThreshold(window, threshold)
    ccall((:sfWindow_setJoystickThreshold, libcsfml_window), Cvoid, (Ptr{sfWindow}, Cfloat), window, threshold)
end

function sfWindow_setActive(window, active)
    ccall((:sfWindow_setActive, libcsfml_window), sfBool, (Ptr{sfWindow}, sfBool), window, active)
end

function sfWindow_requestFocus(window)
    ccall((:sfWindow_requestFocus, libcsfml_window), Cvoid, (Ptr{sfWindow},), window)
end

function sfWindow_hasFocus(window)
    ccall((:sfWindow_hasFocus, libcsfml_window), sfBool, (Ptr{sfWindow},), window)
end

function sfWindow_display(window)
    ccall((:sfWindow_display, libcsfml_window), Cvoid, (Ptr{sfWindow},), window)
end

function sfWindow_getSystemHandle(window)
    ccall((:sfWindow_getSystemHandle, libcsfml_window), sfWindowHandle, (Ptr{sfWindow},), window)
end
# Julia wrapper for header: /Users/gnimuc/.julia/dev/LibCSFML/gen/../deps/usr/include/SFML/Window/WindowHandle.h
# Automatically generated using Clang.jl


text:
@testset "Observations" begin
    @info "Test: Observations"

    @testset "API" begin
        @inferred Observations(tuple(data))
    end
end


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_optimizer/optimize/nonconvex/postprocess.jl
# Routine which calls duality-based bound tightening after solving the upper
# bounding problem.
#############################################################################

"""
$(TYPEDSIGNATURES)

Default postprocess perfoms duality-based bound tightening (Tawarmalani, M.,
Sahinidis, N.V.: Global optimization of mixed-integer nonlinear programs:
a theoretical and computational study. Math. Progr. 99, 563–591 (2004).) up to
an iteration limit set by (`dbbt_depth`).
"""
function postprocess!(t::ExtensionType, m::GlobalOptimizer)
    if m._parameters.dbbt_depth > m._iteration_count
        variable_dbbt!(m._current_node, m._lower_lvd, m._lower_uvd,
                       m._lower_objective_value, m._global_upper_bound,
                       m._branch_variable_count)
    end
    return
end
postprocess!(m::GlobalOptimizer{R,S,Q}) where {R,S,Q<:ExtensionType} = postprocess!(_ext(m), m)

"""
$(TYPEDSIGNATURES)

Checks to see if current node should be reprocessed.
Currently, this returns `false`` by default.
"""
repeat_check(t::ExtensionType, m::GlobalOptimizer) = false
repeat_check(m::GlobalOptimizer{R,S,Q}) where {R,S,Q<:ExtensionType} = repeat_check(_ext(m), m)

max_stars_repo_name: MatthewStuber/EAGO

max_issues_repo_name: MatthewStuber/EAGO

max_forks_repo_name: MatthewStuber/EAGO


text:
# This script prints out HackBio participant data as a part of stage 0 task.
# AUTHOR: Pavlo Hrab

print("NAME: Pavlo Hrab\nE-MAIL: pavlo.hrab@lnu.edu.ua\nUSED LANGUAGE: JULIA\nBIOSTACK: Genomics\nSLACK USERNAME: @pavlo\n")

max_stars_repo_name: jravimicro/team-rosalind-project

max_issues_repo_name: jravimicro/team-rosalind-project

max_forks_repo_name: jravimicro/team-rosalind-project


text:
####################### kernel for SARA's inner loop ###########################
# lx is length scale for position, respectively
function iSARA(lx::Real, rescale_parameters = nothing)
    kx = Kernel.Lengthscale(Kernel.MaternP(2), lx)
    kx = Kernel.SymmetricKernel(kx)
    if !isnothing(rescale_parameters)
        kx = Kernel.VerticalRescaling(kx, get_rescale(rescale_parameters))
    end
    return kx = kx + Kernel.Line(1.)
end
function iSARA(θ::AbstractVector, rescale_parameters = nothing)
    lx = exp(θ[1])
    iSARA(lx, rescale_parameters)
end

############################# rescaling function ###############################
# the scaling function will take the RGB image gradients directly,
# and convolute them with a window function as follows:
# LSA(x) * (grad(RGB(x)^2 * alpha + 1)
# generalized bell curve
gaussian(μ::Real, σ::Real, β::Real) = x->exp(-(abs(x-μ)/σ)^β/2)

# Parameters is a tuple of gaussian parameters, with (c, μ, σ, β)
function sum_gaussians(x::Real, parameters::Vector{<:NTuple{4, T}}) where {T}
    f = zero(x)
    @inbounds @simd for par in parameters
        f += par[1] * gaussian(par[2], par[3], par[4])(x)
    end
    return f
end

function get_rescale(parameters::Vector{<:NTuple{4, <:Real}})
    (x::Real)->sum_gaussians(x, parameters)
end
get_rescale(p::NTuple{4, <:Real}) = get_rescale([p])

################################################################################
# position is from center of stripe
# optical is a vector of vectors of optical coefficients of position
# i.e. optical[1] is the first coefficient as a function of position
# domain is the minimum and maximum position
# rescale are the parameters of the rescaling function
# sigma is the standard deviation (NOT VARIANCE) of the noise in the coefficients
# grid is a pre-determined grid of potential measurement points
# returns: a vector of points to be measured
function inner(position::Vector{<:Real}, optical::Vector{<:Vector},
                domain::NTuple{2, <:Real} = (-1, 1), σ = .01,
                policy = uncertainty_sampling; # sampling policy
                rescale_parameters = nothing,
                grid::AbstractVector{<:Real} = range(domain[1], domain[2], length = 256))
    # data pre-processing
    x = position
    y = optical
    normalization = [(mean(yi), std(yi)) for yi in y]
    y = center.(y) # does not overwrite input data
    # set up prior Gaussian process
    θ = [.2] # length scale, chosen by cross-validation
    log_θ = log.(θ)
    k = iSARA(log.(θ), rescale_parameters)
    G = Gaussian(k)

    # posterior Gaussian processes
    posteriors = Vector{Gaussian}(undef, length(y))
    @threads for i in eachindex(y)
        m, s = normalization[i]
        posteriors[i] = G | (x, y[i], (σ/s)^2)
        posteriors[i] = s*posteriors[i] + m
    end

    # active learning acquisition
    P = posteriors[1]
    next = policy(P, grid)
    next = [next]
    gradient_uncertainty = integrated_uncertainty(P, grid)
    return next, posteriors, gradient_uncertainty
end
const innerloop = inner

########################## local model evaluation ##############################
function evaluate(G::AbstractVector{<:Gaussian}, x::AbstractVector)
    n = length(G)
    T = eltype(x)
    M = Vector{Vector{T}}(undef, n)
    S = Vector{Vector{T}}(undef, n)
    for (i, P) in enumerate(G)
        M[i], S[i] = evaluate(P, x)
    end
    return M, S
end

function evaluate(G::Gaussian, x::AbstractVector)
    M = marginal(G, x)
    return mean(M), std(M)
end


text:
module Parameters

export Param, setParameters

const BCDirichlet,BCRiemann = 1,2 # Constants for different boundary types

using TOML

# Parameters
struct Param
   nD::Int             # Dimension of the system
   LatticeType::Int    # DmQn lattice model
   Δx::Float64         # Discretized space interval
   Δy::Float64         # Discretized space interval (Δy = Δx)
   Δt::Float64         # Discretized time interval
   u ::Vector{Float64} # Advection velocity
   α ::Float64         # Thermal diffusivity
   ω ::Float64         # Relaxation coefficient
   TotalTime::Float64  # Total simulation time
   nT::Int             # Total timesteps
   nI::Int             # Number of space intervals in x
   nJ::Int             # Number of space intervals in y
   T::Vector{Float64}  #  boundary temperatures
   BCtype::Array{Int,1}# Type of boundary conditions
   DoEcho::Bool        # Display parameters on screen

   function Param(nD=2, LatticeType=3, dx=1.0, dy=1.0, dt=1.0, u=[0.1,0.0],
      α=0.25, TotalTime=400.0, nI=100, nJ=100, Tᴸ=1.0, Tᴿ:=0.0, Tᴮ=0.0, Tᵀ=0.0,
      BCtype=[1,1,1,1], DoEcho=true)
      @assert(0 ≤ TotalTime, "Simulation time must be positive!")

      new(nD,LatticeType,dx,dy,dt,u,α,inv(nD*dt*α/dx^2 + 0.5),
         TotalTime,Int(TotalTime/dt),nI,nJ,[Tᴸ,Tᴿ,Tᴮ,Tᵀ],BCtype,DoEcho)
   end
end

function setParameters()
   paramIn = TOML.parsefile("PARAM.toml")
   DoEcho = paramIn["Parameters"]["DoEcho"]
   nD = paramIn["Parameters"]["nD"]
   LatticeType = paramIn["Parameters"]["LatticeType"]
   dx = paramIn["Parameters"]["dx"]
   dy = paramIn["Parameters"]["dy"]
   dt = paramIn["Parameters"]["dt"]
   u = paramIn["Parameters"]["u"]
   α = paramIn["Parameters"]["α"]
   TotalTime = paramIn["Parameters"]["TotalTime"]
   nI = paramIn["Parameters"]["nI"]
   nJ = paramIn["Parameters"]["nJ"]
   Tᴸ = paramIn["Parameters"]["Tleft"]
   Tᴿ = paramIn["Parameters"]["Tright"]
   Tᴮ = paramIn["Parameters"]["Tbot"]
   Tᵀ = paramIn["Parameters"]["Ttop"]

   if LatticeType == "D2Q4"
      LT = 1
   elseif LatticeType == "D2Q5"
      LT = 2
   elseif LatticeType == "D2Q9"
      LT = 3
   end

   nBC = length(paramIn["Parameters"]["BCtype"])
   BCtype = Vector{Int}(undef,nBC)
   for i = 1:nBC
       if paramIn["Parameters"]["BCtype"][i] == "Dirichlet"
           BCtype[i] = 1
       elseif paramIn["Parameters"]["BCtype"][i] == "Riemann"
           BCtype[i] = 2
       end
   end

   param = Param(nD,LT,dx,dy,dt,u,α,TotalTime,nI,nJ,Tᴸ,Tᴿ,Tᴮ,Tᵀ,BCtype,DoEcho)

end

end


text:
import Base.length
export plot_coeffs,plot_coeffs!,plot_eigvals

length(m::LTVStateSpaceModel) = size(m.At,3)
"""
x' = predict(model, d)
x' = predict(model, x, u, t)

Form one-step prediction. If model is an LTVmodel, `x,u` and `model` must have the same length.
If `t` is provided, use model at time `t` to predict a single output only.
"""
function predict(model::LTVStateSpaceModel, d)
    T = length(d)
    n = nstates(d)
    m = ninputs(d)
    x = oftype(Matrix, state(d))
    u = oftype(Matrix, input(d))
    @assert T<=length(model) "Can not predict further than the number of time steps in the model"
    xnew = Array{eltype(x)}(undef,n,T)
    @views for t = 1:T
        xnew[:,t] = (model.At[:,:,t] * x[:,t] + model.Bt[:,:,t] * u[:,t])'
    end
    xnew
end

function predict(model::LTVStateSpaceModel, d::AbstractIdData)
    T = length(d)
    n = nstates(d)
    m = ninputs(d)
    x = oftype(Matrix, state(d))
    @assert T<=length(model) "Can not predict further than the number of time steps in the model"
    xnew = Array{eltype(x)}(undef,n,T)
    @views for t = 1:T
        xnew[:,t] = model.At[:,:,t] * x[:,t]
    end
    xnew
end

function predict(model::LTVStateSpaceModel, x, u, i)
    # @assert i<=length(model) "Can not predict further than the number of time steps in the model"
    xnew = model.At[:,:,i] * x + model.Bt[:,:,i] * u
end

"""
x' = simulate(model, x0, u)

Simulate model forward in time from initial state `x0`. If model is an LTVmodel, `u` and `model` must have the same length.
"""
function simulate(model::LTVStateSpaceModel, x0, u)
    T = length(d)
    n = nstates(d)
    m = ninputs(d)
    u = oftype(Matrix, input(d))
    @assert T > n "The calling convention for u is that time is the second dimention (n,T = size(u))"
    @assert T<=length(model) "Can not simulate further than the number of time steps in the model"
    x = Array{eltype(u)}(undef,n,T)
    x[:,1] = x0
    @views for t = 1:T-1
        x[:,t+1] = model.At[:,:,t] * x[:,t] + model.Bt[:,:,t] * u[:,t]
    end
    x
end

function df(model::LTVStateSpaceModel, x, u)
    fx  = model.At
    fu  = model.Bt
    fxx = []
    fxu = []
    fuu = []
    return fx,fu,fxx,fxu,fuu
end

function plot_eigvals(model::LTVModels.LTVStateSpaceModel, plot_circle=true)
    scatter(hcat([eigvals(model.At[:,:,t]) for t=1:length(model)]...)', title="Eigenvalues of dynamics matrix")
    if plot_circle
        ϕ = range(0, stop=2π, length=100)
        plot!(cos.(ϕ),sin.(ϕ), l=(:black,1))
    end
end

function plot_coeffs!(model::LTVStateSpaceModel; kwargs...)
    n,T = size(model.At,1,3)
    plot!(reshape(model.At,n^2,T)'; title="Coefficients of dynamics matrix", kwargs...)
end

function plot_coeffs(model::LTVStateSpaceModel; kwargs...)
    fig = plot()
    plot_coeffs!(model; kwargs...)
    fig
end


text:
using CUDA: CUDA, @cuda
using DataFrames: DataFrame
using LoopVectorization: @avx

import ProbabilisticCircuits: balance_threads

export ExpectationBit

"""
    ExpectationBit(pc::ProbCircuit, lc::LogisticCircuit, data)
    ExpectationBit(pbc::ParamBitCircuitPair, pc::ProbCircuit, lc::LogisticCircuit, data)

Compute Expected Prediction of a Logistic/Regression Circuit w.r.t to a ProbabilistcCircuit.
Uses the BitCircuitPair implementation, which enables parrallism on both CPU/GPU.
If the data is on gpu the gpu version will get called (repectively for cpu.)

Missing values should be denoted by missing
See: On Tractable Computation of Expected Predictions [arxiv.org/abs/1910.02182](https://arxiv.org/abs/1910.02182)
"""
function ExpectationBit(pc::ProbCircuit, lc::LogisticCircuit, data; return_aux = false)
    # children(lc)[1]: Currently root of LC's are for bias values and not a real sum node, so should be skipped.
    pbc = ParamBitCircuitPair(pc, children(lc)[1]);
    pbc = same_device(pbc, data)
    ExpectationBit(pbc, pc, lc, data; return_aux);
end

function ExpectationBit(pbc::ParamBitCircuitPair, pc::ProbCircuit, lc::LogisticCircuit, data, reuse_f=nothing, reuse_g=nothing; return_aux=true)
    # 1. Get probability of each observation
    parambit::ParamBitCircuit = ParamBitCircuit(pbc.pc_bit, pbc.pc_params);
    log_likelihoods = MAR(parambit, data);
    p_observed = exp.( log_likelihoods )

    # 2. Expectation w.r.t. P(x_m, x_o)
    fvalues, gvalues = init_expectations(pbc, data, reuse_f, reuse_g, size(pbc.pair_bit.nodes)[2], num_classes(lc))
    expectation_layers(pbc, fvalues, gvalues)

    # 3. Expectation w.r.t P(x_m | x_o)
    results = gvalues[:, end,:] ./ p_observed

    # # 4. Add Bias terms
    biases = if isgpu(results)
        biases = to_gpu(lc.thetas)
    else
        biases = lc.thetas
    end
    results .+= biases

    if return_aux
        results, fvalues, gvalues, pbc
    else
        results
    end
end

"""
    init_expectations(circuit::ParamBitCircuitPair, data, reuse_f, reuse_g, nodes_num, classes_num; Float = Float32)

    Initialized the input layer for `fvalues` and `gvalues`.
    To reduce memory allocation can pass `reuse_f` and `reuse_g` to get reused (very useful during batching).
"""
function init_expectations(circuit::ParamBitCircuitPair, data, reuse_f, reuse_g, nodes_num, classes_num; Float = Float32)
    flowtype = isgpu(data) ? CuArray{Float} : Array{Float}
    @assert isgpu(data) == isgpu(circuit)

    fvalues = similar!(reuse_f, flowtype, num_examples(data), nodes_num)
    fgvalues = similar!(reuse_g, flowtype, num_examples(data), nodes_num, classes_num)
    nfeatures = num_features(data)

    # TODO might not need to zero everything out, but was giving wrong answer
    # on gpu reuse if not;
    fvalues[:,:] .= zero(Float)
    fgvalues[:,:,:] .= zero(Float)

    PARS = circuit.lc_params
    for var=1:nfeatures
        # even if data is on gpu this might become on cpu (eg when data is dataframe)
        temp_cpu = .!isequal.(data[:, var], 0) .* one(Float)
        fvalues[:, var] .= same_device(temp_cpu, data)

        temp_cpu .= .!isequal.(data[:, var], 1) .* one(Float)
        fvalues[:, var + 3*nfeatures] .= same_device(temp_cpu, data)
    end

    if isgpu(circuit)
        examples_count = size(fvalues, 1)
        num_threads = balance_threads(examples_count, nfeatures/8, 8)
        num_blocks = (ceil(Int, examples_count/num_threads[1]),
                        ceil(Int, nfeatures/num_threads[2]))
        # CUDA.@sync begin
        #     @cuda threads=num_threads blocks=num_blocks init_expectations_cuda_part1(nfeatures, data, fvalues)
        # end

        CUDA.@sync begin
            @cuda threads=num_threads blocks=num_blocks init_expectations_cuda_part2(nfeatures, classes_num,
                                                        circuit.lc_bit.nodes, circuit.lc_bit.parents,
                                                        PARS, fvalues, fgvalues)
        end
    else
        for var=1:nfeatures
            ind2 = var + 3*nfeatures
            # find the index of correct paramter to grab from LogisticCircuit
            p_ind = circuit.lc_bit.parents[circuit.lc_bit.nodes[3, 2+var]]
            p_ind2 = circuit.lc_bit.parents[circuit.lc_bit.nodes[3, 2+var+nfeatures]]

            for cc=1:classes_num
                fgvalues[:, var, cc] .= (fvalues[:, var] .* PARS[p_ind, cc])
                fgvalues[:, ind2, cc] .= (fvalues[:, ind2] .* PARS[p_ind2, cc])
            end
        end
    end
    return fvalues, fgvalues
end

function init_expectations_cuda_part2(nfeatures, classes_num, lc_bit_nodes, lc_bit_pars, PARS, fvalues, fgvalues)
    index_x = (blockIdx().x - 1) * blockDim().x + threadIdx().x
    index_y = (blockIdx().y - 1) * blockDim().y + threadIdx().y
    stride_x = blockDim().x * gridDim().x
    stride_y = blockDim().y * gridDim().y

    for var = index_x:stride_x:nfeatures
        @inbounds ind2 = var + 3*nfeatures
        @inbounds idx1 = lc_bit_nodes[3, 2+var]
        @inbounds idx2 = lc_bit_nodes[3, 2+var+nfeatures]

        @inbounds p_ind = lc_bit_pars[idx1]
        @inbounds p_ind2 = lc_bit_pars[idx2]
        for i = index_y:stride_y:size(fvalues, 1)
            for cc=1:classes_num
                fgvalues[i, var, cc] = (fvalues[i, var] * PARS[p_ind, cc])
                fgvalues[i, ind2, cc] = (fvalues[i, ind2] * PARS[p_ind2, cc])
            end
        end
    end
end

function expectation_layers(circuit, fvalues::Array{<:AbstractFloat,2}, fgvalues::Array{<:AbstractFloat,3})
    bcp::BitCircuitPair = circuit.pair_bit
    pc::BitCircuit = circuit.pc_bit
    lc::BitCircuit = circuit.lc_bit
    els = bcp.elements
    pc_pars = exp.(circuit.pc_params)
    lc_pars = circuit.lc_params

    for layer in bcp.layers[2:end]
        Threads.@threads for dec_id in layer

            id_begin = @inbounds bcp.nodes[1, dec_id]
            id_end   = @inbounds bcp.nodes[2, dec_id]
            pc_id    = @inbounds bcp.nodes[5, dec_id]
            lc_id    = @inbounds bcp.nodes[6, dec_id]
            # println("!!, $dec_id,  $id_begin, $id_end, $pc_id, $lc_id")

            @inbounds pc_childs = 1 + pc.nodes[2, pc_id] - pc.nodes[1, pc_id]
            @inbounds lc_childs = 1 + lc.nodes[2, lc_id] - lc.nodes[1, lc_id]
            for A =  1 : pc_childs
                for B = 1 : lc_childs
                    shift = (A-1) * lc_childs + B
                    Z = id_begin + shift - 1
                    @inbounds PCZ = pc.nodes[1, pc_id] + A - 1
                    @inbounds LCZ = lc.nodes[1, lc_id] + B - 1

                    ## Compute fvalues[:, dec_id]
                    # @inbounds @fastmath @avx fvalues[:, dec_id] .+= pc_pars[PCZ] .*
                    #     fvalues[:, els[2,Z]] .*
                    #     fvalues[:, els[3,Z]]

                    # The for loops are faster, cause much less memory needed
                    # @avx throws error, not sure why
                    for j=1:size(fvalues,1)
                        @inbounds fvalues[j, dec_id] += pc_pars[PCZ] * fvalues[j, els[2,Z]] * fvalues[j, els[3,Z]]
                    end

                    ## Compute fgvalues[:, dec_id, :]
                    if els[3,Z] == els[2,Z]
                        # Special case (end in sum nodes)
                        @inbounds @fastmath @views @avx fgvalues[:, dec_id, :] .=
                            transpose(lc_pars[LCZ,:]) .*
                            fvalues[:, dec_id]
                    else
                        # @inbounds @fastmath @avx fgvalues[:, dec_id, :] .+=
                        #     pc_pars[PCZ] .*
                        #     transpose(lc_pars[LCZ,:]) .*
                        #     (fvalues[:, els[2,Z]] .* fvalues[:, els[3,Z]])

                        # @inbounds @fastmath @avx fgvalues[:, dec_id, :] .+=
                        #     pc_pars[PCZ] .*
                        #     ((fvalues[:, els[3,Z]] .* fgvalues[:, els[2,Z], :]) .+
                        #     (fvalues[:, els[2,Z]] .* fgvalues[:, els[3,Z], :]))

                        # The for loops are faster, cause much less memory needed
                        # @avx throws error, not sure why
                        @simd for j=1:size(fgvalues,1)
                            @simd for nc=1:size(fgvalues, 3)
                                @inbounds @fastmath fgvalues[j, dec_id, nc] +=
                                    pc_pars[PCZ] *
                                    lc_pars[LCZ,nc] *
                                    (fvalues[j, els[2,Z]] * fvalues[j, els[3,Z]])

                                @inbounds @fastmath fgvalues[j, dec_id, nc] +=
                                    pc_pars[PCZ] *
                                    ((fvalues[j, els[3,Z]] * fgvalues[j, els[2,Z], nc]) +
                                    (fvalues[j, els[2,Z]] * fgvalues[j, els[3,Z], nc]))
                            end
                        end
                    end

                end # B
            end # A

        end # dec_id
    end # layers

end

function expectation_layers(circuit::ParamBitCircuitPair,
    fvalues::CuArray, fgvalues::CuArray;
    dec_per_thread = 8, log2_threads_per_block = 8)

    bcp = circuit.pair_bit
    pc_pars = exp.(circuit.pc_params)

    CUDA.@sync for layer in bcp.layers[2:end]
        num_examples = size(fvalues, 1)
        num_decision_sets = length(layer)/dec_per_thread
        num_threads =  balance_threads(num_examples, num_decision_sets, log2_threads_per_block)
        num_blocks = (ceil(Int, num_examples/num_threads[1]),
                      ceil(Int, num_decision_sets/num_threads[2]))

        @cuda threads=num_threads blocks=num_blocks expectation_layers_cuda(layer, bcp.nodes, bcp.elements,
            circuit.pc_bit.nodes, circuit.lc_bit.nodes,
            pc_pars, circuit.lc_params, fvalues, fgvalues)
    end
    return nothing
end

function expectation_layers_cuda(layer, nodes, els,
    pc_nodes, lc_nodes,
    pc_pars, lc_pars, fvalues, fgvalues)

    num_classes = size(lc_pars, 2)

    index_x = (blockIdx().x - 1) * blockDim().x + threadIdx().x
    index_y = (blockIdx().y - 1) * blockDim().y + threadIdx().y
    stride_x = blockDim().x * gridDim().x
    stride_y = blockDim().y * gridDim().y
    for j = index_x:stride_x:size(fvalues,1)
        for i = index_y:stride_y:length(layer)
            @inbounds dec_id =  layer[i]
            @inbounds id_begin  =  nodes[1, dec_id]
            @inbounds id_end    =  nodes[2, dec_id]
            @inbounds pc_id     =  nodes[5, dec_id]
            @inbounds lc_id     =  nodes[6, dec_id]
            @inbounds pc_childs = 1 + pc_nodes[2, pc_id] - pc_nodes[1, pc_id]
            @inbounds lc_childs = 1 + lc_nodes[2, lc_id] - lc_nodes[1, lc_id]

            for A =  1 : pc_childs
                for B = 1 : lc_childs
                    shift = (A-1) * lc_childs + B
                    Z = id_begin + shift - 1
                    @inbounds PCZ = pc_nodes[1, pc_id] + A - 1
                    @inbounds LCZ = lc_nodes[1, lc_id] + B - 1

                    ## Compute fvalues[:, dec_id]
                    fvalues[j, dec_id] += pc_pars[PCZ] *
                        fvalues[j, els[2,Z]] *
                        fvalues[j, els[3,Z]]

                    ## Compute fgvalues[:, dec_id, :]
                    if els[3,Z] == els[2,Z]
                        # Special case (end in sum nodes)
                        for class=1:num_classes
                            @inbounds fgvalues[j, dec_id, class] =
                                lc_pars[LCZ,class] *
                                fvalues[j, dec_id]
                        end
                    else
                        for class=1:num_classes
                            @inbounds fgvalues[j, dec_id, class] +=
                                pc_pars[PCZ] *
                                lc_pars[LCZ,class] *
                                (fvalues[j, els[2,Z]] * fvalues[j, els[3,Z]])

                            @inbounds fgvalues[j, dec_id, class] +=
                                pc_pars[PCZ] *
                                ((fvalues[j, els[3,Z]] * fgvalues[j, els[2,Z], class]) +
                                (fvalues[j, els[2,Z]] * fgvalues[j, els[3,Z], class]))

                        end
                    end
                end # B
            end # A
        end # i
    end # j
    return nothing
end

max_stars_repo_name: Juice-jl/DiscriminativeCircuits

max_issues_repo_name: Juice-jl/DiscriminativeCircuits

max_forks_repo_name: Juice-jl/DiscriminativeCircuits


text:
# ## More data processing
#
# This tutorial uses the World Resources Institute Global Power Plants Dataset to explore data pre-processing in Julia.
# The dataset is created from multiple sources and is under continuous update, which means that there are lots of missing data, non-standard characters, etc
# Hence plenty of material to work with!

# More tutorials on the manipulation of DataFrames can be found [here](https://github.com/bkamins/Julia-DataFrames-Tutorial)
# And some more information can be found on [this](https://en.wikibooks.org/wiki/Introducing_Julia/DataFrames) wikipage.

import MLJ: schema, std, mean, median, coerce, coerce!, scitype
using DataFrames
using UrlDownload
using PyPlot
ioff() # hide

# Import data

raw_data = urldownload("https://github.com/tlienart/DataScienceTutorialsData.jl/blob/master/data/wri_global_power_plant_db_be_022020.csv?raw=true")
data = DataFrame(raw_data);

# This dataset contains information on power generation plants for a number of countries around the world.
# The level of disaggregation is the power plant. For each plant, there is information about its name, localisation, capacity, and many other features.
# The schema function enables us to get a quick overview of the variables it contains, including their machine and scentific types.

schema(data)

# We see that a small number of features have values for all plants (i.e. for each row) present in the dataset.
# However, (i) several features have missing values (Union{Missing, _.type}) and (ii) we are not interested in working with all of these features.
# In particular, we're not intersted in the source of the information present in the dataset nor are we interested in the generation data.
# Hence we drop all columns which contain information's source.
# We define a function `is_active()` that will return a `TRUE` boolean value if the column name does NOT (`!`) contain either of the strings "source" or "generation".
# Note the conversion of column names from `:Symbol` to `:string` since the `occursing` function only accepts strings as arguments.

is_active(col) = !occursin(r"source|generation", string(col))
active_cols = [col for col in names(data) if is_active(col)]
select!(data, active_cols);

# We also drop a number of other unwanted columns and take a look at our "new" dataframe.

select!(data, Not([:wepp_id, :url, :owner]))
schema(data)

# The remaining variables have two different scientific types: Continuous, Textual
# Of which we can get an overview.

describe(data)

# The describe() function shows that there are several features with missing values.
# *Note:* the `describe()` function is from the `DataFrames` package (and won't work with other, non DataFrames, tables) whereas the `schema()` is from the MLJ package.

# ---
# Let's play around with capacity data, for which there are no missing values. We create a sub-dataframe and aggregate over certain dimensions (country and primary_fuel)
capacity = select(data, [:country, :primary_fuel, :capacity_mw]);
first(capacity, 5)

# This dataframe contains several subgroups (country and technology type) and it would be interesting to get data aggregates by subgroup.
# To obtain a `view` of the DataFrame by subgroup, we can use the `groupby` function.
# (See the [DataFrame tutorial](https://alan-turing-institute.github.io/DataScienceTutorials.jl/data/dataframe/#groupby) for an introduction to the use of `groupby`)
cap_gr = groupby(capacity, [:country, :primary_fuel]);

# If we want to aggregate at the country-fuel-type level and calculate summary statistics at this level, we can use the `combine` function on the GroupedDataFrame that we just created.
# This function takes the GroupedDataFrame, the symbol of the column on which to apply the measure of choice as arguments.

cap_mean = combine(cap_gr, :capacity_mw => mean)
cap_sum = combine(cap_gr, :capacity_mw => sum)
first(cap_sum, 3)

# Now let's plot some of this aggregate data for a selection of countries, by country and technology type

ctry_selec = r"BEL|FRA|DEU"
tech_selec = r"Solar"

cap_sum_plot = cap_sum[occursin.(ctry_selec, cap_sum.country) .& occursin.(tech_selec, cap_sum.primary_fuel), :]

# Note the `.` for element-wise comparison
# Before plotting, we can also sort values by decreasing order using `sort!()`.
sort!(cap_sum_plot, :capacity_mw_sum, rev=true)

figure(figsize=(8,6))

plt.bar(cap_sum_plot.country, cap_sum_plot.capacity_mw_sum, width=0.35)
plt.xticks(rotation=90)

savefig(joinpath(@OUTPUT, "D0-processing-g1.svg")) # hide

# \figalt{processing1}{D0-processing-g1.svg}

# ---
# Now that we have the total capacity by country and technology type, let's use it to calculate the share of each technology in total capacity.
# To that end we first create a dataframe containing the country-level total capacity, using the same steps as above.
cap_sum_ctry_gd = groupby(capacity, [:country]);
cap_sum_ctry = combine(cap_sum_ctry_gd, :capacity_mw => sum);

# The we join this dataframe with the disaggregated one; which requires that we convert the two GroupedDataFrame into DataFrames.

cap_sum = DataFrame(cap_sum);
cap_sum_ctry = DataFrame(cap_sum_ctry);
cap_share = leftjoin(cap_sum, cap_sum_ctry, on = :country, makeunique = true)
cap_share.capacity_mw_share = cap_share.capacity_mw_sum ./ cap_share.capacity_mw_sum_1;

# Let's visualise our dataframe again, which now includes the `capacity_mw_share` column.

# ---
# Now let's analyse features which exhibit some missing values.
# Suppose we want to calculate the age of each plant (rounded to full years). We face two issues.
# First, the commissioning_year is not reported for all plants.
# We need to gauge the representativity of the plants for which it is available with regard to the full dataset.
# One way to count the missing values is

nMissings = length(findall(x -> ismissing(x), data.commissioning_year))

# This represents about half of our observations
nMissings_share = nMissings/size(data)[1]

# Second, the commissioning year is not reported as an integer. Fractions of years are reported too.
# As a result, the machine type of `data.commissioning_year`is Float64.

typeof(data.commissioning_year)

# Before calculating the average age, let's drop the missing values.

data_nmiss = dropmissing(data, :commissioning_year);

# And round the year to the closest integer. We can do this using the `round` function and a mapping function on the relevant DataFrame column.

map!(x -> round(x, digits=0), data_nmiss.commissioning_year, data_nmiss.commissioning_year);

# We can now calculate plant age for each plant (worth remembering that the dataset only contains active plants)

current_year = fill!(Array{Float64}(undef, size(data_nmiss)[1]), 2020);
data_nmiss[:, :plant_age] = current_year - data_nmiss[:, :commissioning_year];

# Since the commissioning year is missing for about half the plants in the dataset (17340, see description of data above) and that missing values propagate,
# the plant age will only be available for 33643-17340 plants.
# Let's see what the mean and median plant ages are across the plants for which we have the data

mean_age = mean(skipmissing(data_nmiss.plant_age))
median_age = median(skipmissing(data_nmiss.plant_age))

# And bring this into a frequency plot of the plant age observations

figure(figsize=(8,6))

plt.hist(data_nmiss.plant_age, color="blue", edgecolor="white", bins=100,
      density=true, alpha=0.5)
plt.axvline(mean_age, label = "Mean", color = "red")
plt.axvline(median_age, label = "Median")

plt.legend()

plt.xlim(0,)

savefig(joinpath(@OUTPUT, "D0-processing-g2.svg")) # hide

# \figalt{processing2}{D0-processing-g2.svg}

# We can also calculate and plot average plant age by country and technology
# Make sure all columns passed, other than the aggregation dimensions, are of type `Float` or `Int`, otherwise the function execution will fail.

age = select(data_nmiss, [:country, :primary_fuel, :plant_age])
age_mean = combine(groupby(age, [:country, :primary_fuel]), :plant_age => mean)

coal_means = age_mean[occursin.(ctry_selec, age_mean.country) .& occursin.(r"Coal", age_mean.primary_fuel), :]
gas_means = age_mean[occursin.(ctry_selec, age_mean.country) .& occursin.(r"Gas", age_mean.primary_fuel), :]

width = 0.35  # the width of the bars

fig, (ax1, ax2) = plt.subplots(1,2)

fig.suptitle("Mean plant age by country and technology")

ax1.bar(coal_means.country, coal_means.plant_age_mean, width, label="Coal")
ax2.bar(gas_means.country, gas_means.plant_age_mean, width, label="Gas")

ax1.set_ylabel("Age")

ax1.set_title("Coal")
ax2.set_title("Gas")

savefig(joinpath(@OUTPUT, "D0-processing-g3.svg")) # hide

# \figalt{processing3}{D0-processing-g3.svg}


text:
function khayyam_triangle(lentingle)
    output = []
    for i in 1:lentingle
        push!(output, [])
        for j in 1:i
            if j == 1
                append!(output[i], 1)
            elseif j == i
                append!(output[i], 1)
            else j != 1 && j != i
                add = output[i - 1][j] + output[i - 1][j - 1]
                append!(output[i], add)
            end
        end
    end
    return output
end

max_stars_repo_name: shahmari/A-few-fractals-in-Julia

max_issues_repo_name: shahmari/A-few-fractals-in-Julia

max_forks_repo_name: shahmari/A-few-fractals-in-Julia


text:
@with_kw mutable struct UrbanDriver <: DriverModel{LonAccelDirection}
    a::LonAccelDirection = LonAccelDirection(0., 1)
    navigator::RouteFollowingIDM = RouteFollowingIDM()
    intersection_driver::DriverModel = StopIntersectionDriver()
    crosswalk_drivers::Vector{CrosswalkDriver} = [CrosswalkDriver()]

    debug::Bool = false
end

function AutomotiveDrivingModels.reset_hidden_state!(model::UrbanDriver)
    reset_hidden_state!(model.navigator)
    reset_hidden_state!(model.intersection_driver)
    for driver in model.crosswalk_drivers
        reset_hidden_state!(driver)
    end
    model
end

function AutomotiveDrivingModels.observe!(model::UrbanDriver, scene::Scene, roadway::Roadway, egoid::Int)
    AutomotiveDrivingModels.observe!(model.navigator, scene, roadway, egoid)
    AutomotiveDrivingModels.observe!(model.intersection_driver, scene, roadway, egoid)
    for driver in model.crosswalk_drivers
        AutomotiveDrivingModels.observe!(driver, scene, roadway, egoid)
    end
    a_lon_crosswalks = minimum([driver.a.a_lon for driver in model.crosswalk_drivers])
    # println("intersection driver ", model.intersection_driver.a)
    # println("crosswalk driver ", model.crosswalk_driver.a)
    # take the more conservative of the three actions
    a_lon = min(model.intersection_driver.a.a_lon, a_lon_crosswalks)
    # a_lon = model.intersection_driver.a.a_lon
    # i = indmin([model.intersection_driver.a.a_lon, model.crosswalk_drivers[1].a.a_lon, model.crosswalk_drivers[2].a.a_lon, model.crosswalk_drivers[3].a.a_lon])
    # println(egoid, " getting minimum acceleration from ", i)
    if model.debug
        println("inter ", model.intersection_driver.a.a_lon)
        println("crosswalk 1 ", model.crosswalk_drivers[1].a.a_lon)
    end
    model.a = LonAccelDirection(a_lon, model.navigator.dir)
    return model
end

get_name(::UrbanDriver) = "UrbanDriver"
Base.rand(model::UrbanDriver) = model.a

function action_space(model::UrbanDriver)
    return (LonAccelDirection(model.a.a_lon-model.navigator.σ, model.navigator.dir), LonAccelDirection(model.a.a_lon, model.navigator.dir), LonAccelDirection(model.a.a_lon + model.navigator.σ, model.navigator.dir))
end

function get_distribution(model::UrbanDriver)
    actions = action_space(model)
    probs = (1/3, 1/3, 1/3)
    return actions, probs
end

max_stars_repo_name: MaximeBouton/AutomotivePOMDPs

max_issues_repo_name: MaximeBouton/AutomotivePOMDPs

max_forks_repo_name: MaximeBouton/AutomotivePOMDPs


text:
function test_markov_chain()

	stpm = [0.7 0.2 0.1; 0.3 0.5 0.2; 0 0 1]

	MarkovChainFamily(stpm)

end


text:
# ---
# title: 745. Prefix and Suffix Search
# id: problem745
# author: Tian Jun
# date: 2020-10-31
# difficulty: Hard
# categories: Trie
# link: <https://leetcode.com/problems/prefix-and-suffix-search/description/>
# hidden: true
# ---
#
# Design a special dictionary which has some words and allows you to search the
# words in it by a prefix and a suffix.
#
# Implement the `WordFilter` class:
#
#   * `WordFilter(string[] words)` Initializes the object with the `words` in the dictionary.
#   * `f(string prefix, string suffix)` Returns _the index of the word in the dictionary_ which has the prefix `prefix` and the suffix `suffix`. If there is more than one valid index, return **the largest** of them. If there is no such word in the dictionary, return `-1`.
#
#
#
# **Example 1:**
#
#
#
#     **Input**
#     ["WordFilter", "f"]
#     [[["apple"]], ["a", "e"]]
#     **Output**
#     [null, 0]
#
#     **Explanation**
#     WordFilter wordFilter = new WordFilter(["apple"]);
#     wordFilter.f("a", "e"); // return 0, because the word at index 0 has prefix = "a" and suffix = 'e".
#
#
#
#
# **Constraints:**
#
#   * `1 <= words.length <= 15000`
#   * `1 <= words[i].length <= 10`
#   * `1 <= prefix.length, suffix.length <= 10`
#   * `words[i]`, `prefix` and `suffix` consist of lower-case English letters only.
#   * At most `15000` calls will be made to the function `f`.
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end


text:
"""
    path(multidoc, binding, sig) -> String
Returns the path of this binding with given signature.
"""
path(m::MultiDoc, b::Binding, sig) = m.docs[sig].data[:module] => m.docs[sig].data[:path]

"""
    paths(multidoc, binding) -> Generator
Returns an iterator of path of the documents of this binding.
"""
paths(m::MultiDoc, b::Binding) = (path(m, b, sig) for sig in m.order)

"""
    paths(module) -> Set
Returns an `Set` object that contains the unique relative path of doc strings.
"""
function paths(__module__::Module)
    META = Docs.meta(__module__)
    PATH = Set()
    for (b, m) in META
        for each in paths(m, b)
            push!(PATH, each)
        end
    end
    PATH
end

"""
    dump_str(binding[, sig]) -> String
dump a `binding` and a signature `sig` to raw script.
"""
function dump_str end

dump_str(b::Binding) = "$(b.var)"

function dump_str(b::Binding, sig)
    if sig === Union{}
        "$(b.var)"
    else
        "$(b.var) $sig"
    end
end

function dump_str(b::Binding, sig::Type{<:Tuple})
    o = "$(b.var)"

    args = join(["::$each" for each in sig.types], ", ")
    !isempty(args) && return o * "(" * args * ")"
    o
end

function allpaths(__module__::Module)
    PATH = paths(__module__)
    for each in submodules(__module__)
        union!(PATH, allpaths(each))
    end
    PATH
end

template(mod::Module; content="# translation goes here") =
###### Template
"""
@i18n $mod begin
$content
end
"""

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

readme() = """
# I18n generated files
Files and dirs in this folder is generated by [Localize.jl](). This keeps the original directory structure and adds
i18n macros to each file. Just simply put your translation in each generated Julia script.
"""

include_string(path::String) = "include(\"$path\")"

function generate(__module__::Module, root; workdir=pwd())
    if ispath(root)
        @warn "overwriting existing path: $root"
    else
        mkpath(root)
    end

    open(joinpath(root, "README.md"), "w+") do f
        write(f, readme())
    end

    include_paths = String[]
    files = Dict{String, Vector{Module}}()

    # merge different module paths
    for (mod, each) in allpaths(__module__)
        if haskey(files, each)
            push!(files[each], mod)
        else
            files[each] = [mod]
        end
    end

    for (each, mods) in files
        path = joinpath(root, each)
        push!(include_paths, each)
        mkpath(dirname(path))

        printstyled("path: "; bold=true, color=:light_black)
        println(path)

        open(path, "a+") do f
            printstyled("module: "; bold=true, color=:light_blue)
            for mod in mods
                printstyled(mod; bold=true, color=:yellow)
                print(" ") # make a new line

                write(f, template(mod))
            end
        end
        println("\n")
    end

    include_script = join(include_string.(include_paths), "\n")

    open(joinpath(root, "__main__.jl"), "w+") do f
        write(f, include_script)
        write(f, "\n") # add line break at the EOF
    end
    nothing
end

generate(root) = generate(Base, root)


text:
"""
    WhittleLikelihood(model::Type{<:TimeSeriesModel}, ts, Δ; lowerΩcutoff, upperΩcutoff, taper)
    WhittleLikelihood(model::Type{<:TimeSeriesModel}, timeseries::TimeSeries; lowerΩcutoff, upperΩcutoff, taper)

Generate a function to evaluate the Whittle likelihood it's gradient and Hessian.

Create a callable struct which prealloctes memory appropriately.

    (f::WhittleLikelihood)(θ)

Evaluates the Whittle likelihood at θ.

    (f::WhittleLikelihood)(F,G,H,θ)

Evaluates the Whittle likelihood at θ and stores the gradient and Hessian in G and H respectively.
If F, G or H equal nothing, then the function, gradient or Hessian are not evaluated repsectively.

# Aruguments
- `model`: the model for the process. Should be of type TimeSeriesModel, so `OU` and not `OU(1,1)`.
- `ts`: the timeseries in the form of an n by d matrix (where d is the dimension of the time series model).
- `Δ`: the sampling rate of the time series.
- `timeseries`: can be provided instead of `ts` and `Δ`. Must be of type `TimeSeries`.
- `lowerΩcutoff`: the lower bound of the frequency range included in the likelihood.
- `upperΩcutoff`: the upper bound of the frequency range included in the likelihood.
- `taper`: optional taper which should be a vector of length `size(ts,1)`, or `nothing` in which case no taper will be used.

Note that to use the gradient the model must have `grad_add_sdf!` specified.
Similarly, to use the Hessian, the model must have `hess_add_sdf!` specified.

# Examples
```julia-repl
julia> obj = WhittleLikelihood(OU, ones(1000), 1)
Whittle likelihood for the OU model.

julia> obj([1.0, 1.0])
-2006.7870804551364

julia> F, G, H = 0.0, zeros(2), zeros(2,2)
(0.0, [0.0, 0.0], [0.0 0.0; 0.0 0.0])

julia> obj(F, G, H, [1.0, 1.0])
-2006.7870804551364

julia> G
2-element Vector{Float64}:
   2.777354965282642
 -17.45063591068618

julia> H
2×2 Matrix{Float64}:
 -0.00967179   0.0607696
  0.0607696   -0.381827
```

"""
struct WhittleLikelihood{T,S<:TimeSeriesModelStorage,M}
    data::WhittleData{T}
    memory::S
    function WhittleLikelihood(model, timeseries::TimeSeries; lowerΩcutoff = 0, upperΩcutoff = Inf, taper = nothing)
        WhittleLikelihood(model, timeseries.ts, timeseries.Δ; lowerΩcutoff = lowerΩcutoff, upperΩcutoff = upperΩcutoff, taper = taper)
    end
    function WhittleLikelihood(
        model::Type{<:TimeSeriesModel{D}}, ts, Δ;
        lowerΩcutoff = 0, upperΩcutoff = Inf, taper = nothing) where {D}

        Δ > 0 || throw(ArgumentError("Δ should be a positive."))
        D == size(ts,2) || throw(ArgumentError("timeseries is $(size(ts,2)) dimensional, but model is $D dimensional."))
        if taper !== nothing
            taper .*= inv(sqrt(sum(abs2,taper))) # normalise the taper
        end
        wdata = WhittleData(model, ts, Δ, lowerΩcutoff = lowerΩcutoff, upperΩcutoff = upperΩcutoff, taper = taper)
        mem = allocate_memory_sdf_FGH(model, size(ts,1), Δ)
        new{eltype(wdata.I),typeof(mem),model}(wdata,mem)
    end
end
(f::WhittleLikelihood{T,S,M})(θ) where {T,S,M} = whittle!(f.memory,M(θ),f.data)
function checksize(G,H,θ)
    G === nothing || size(G) == size(θ) || throw(ArgumentError("G should be nothing or same size as θ"))
    H === nothing || size(H) == (length(θ),length(θ)) || throw(ArgumentError("H should be nothing or size (length(θ),length(θ))"))
    nothing
end
function (f::WhittleLikelihood{T,S,M})(F,G,H,θ) where {T,S,M}
    checksize(G,H,θ)
    whittle_FGH!(F,G,H,f.memory,M(θ),f.data)
end
Base.show(io::IO, ::WhittleLikelihood{T,S,M}) where {T,S,M} = print(io, "Whittle likelihood for the $M model.")
getmodel(::WhittleLikelihood{T,S,M}) where {T,S,M} = M
## internal functions

"""
    whittle!(store, model::TimeSeriesModel, data::GenWhittleData)

Compute the Whittle likelihood using a preallocated store.
"""
function whittle!(store, model::TimeSeriesModel, data::GenWhittleData)
    asdf!(store, model)
    return @views generalwhittle(store, data)
end

"""
    whittle_FG!(F, G, store, model::TimeSeriesModel, data::GenWhittleData)

Compute the Whittle likelihood and its gradient using a preallocated store.
"""
function whittle_FG!(F, G, store, model::TimeSeriesModel, data::GenWhittleData)
    if F !== nothing || G !== nothing
        asdf!(model, store)
    end
    if G !== nothing
        grad_asdf!(store, model)
        grad_generalwhittle!(G, store, data)
    end
    if F !== nothing
        return generalwhittle(store, data)
    end
    return nothing
end

"""
    whittle_FGH!(F, G, H, store, model::TimeSeriesModel, data::GenWhittleData)

Compute the debiased Whittle likelihood and its gradient and hessian using a preallocated store.
"""
function whittle_FGH!(F, G, H, store, model::TimeSeriesModel, data::GenWhittleData)
    if F !== nothing || G !== nothing || H !== nothing
        asdf!(store, model)
    end
    if G !== nothing || H !== nothing
        grad_asdf!(store, model)
    end
    if H !== nothing
        hess_asdf!(store, model)
        hess_generalwhittle!(H, store, data)
    end
    if G !== nothing
        grad_generalwhittle!(G, store, data)
    end
    if F !== nothing
        return generalwhittle(store, data)
    end
    return nothing
end


text:
# This script includes a demo for the errot plots.

using Makie
using FractalTools
using GeometryBasics
using Random

Random.seed!(0)

# Generate data
trig = Triangle(
    Point(BigFloat(-1.), BigFloat(-1.)),
    Point(BigFloat(1.), BigFloat(-1)),
    Point(BigFloat(0.), BigFloat(2.))
)
f(x, y) = x^2 + y^2
pnts3d = getdata(f, trig, 100)

# Construct interpolant
itp = interpolate(pnts3d)
function erf(x, y)
    ival = itp(x, y)
    fval = f(x, y)
    abs(fval - ival) / abs(fval) * 100
end

# Construct gui
_, msh3 = triangulate(pnts3d)
msh2 = project(msh3)
pnts2d = project(pnts3d)

fig         = Figure()
lb          = Label(fig, "Show Interpolant")
tog         = Toggle(fig, active=false)
ax1         = LScene(fig, scenekw=(camema=cam3d!, raw=false), tellwidth=false)
fig[1, 1]   = vgrid!(ax1, hgrid!(lb, tog, tellwidth=false), tellwidth=false)
ax2         = fig[1, 2] = LScene(fig, scenekw=(camema=cam3d!, raw=false), tellwidth=false)

        trisurf!(ax1, msh2, f, meshcolor3=first.(msh3.position), colormap=:viridis)
plt =   trisurf!(ax1, msh2, itp, meshcolor3=first.(msh3.position), colormap=:heat)
        wireframe!(ax1, msh2)
        trisurf!(ax2, msh2, erf, meshcolor3=first.(msh2.position))

on(tog.active) do val
    plt.visible[] = val
end

tog.active[] = false

display(fig)

# fps = 60
# record(fig.scene, "errorplot.mp4"; framerate = fps) do io
#     for i = 1 : 10 * fps
#         sleep(1/fps)
#         recordframe!(io)
#     end
# end


text:
using Powseeker
using POMDPToolbox
using POMCP
using POMCPOW
using Plots
using Parameters
using POMDPs
using ContinuousPOMDPTreeSearchExperiments
using ParticleFilters
using Reel
using ProgressMeter

# using Gallium
# Gallium.breakpoint(Pkg.dir("ContinuousPOMDPTreeSearchExperiments","src","ContinuousPOMDPTreeSearchExperiments.jl"), 68)
# @conditional Gallium.breakpoint(Pkg.dir("POMCPOW","src","solver2.jl"), 112) (h_node.node==14)

pomdp = PowseekerPOMDP()
rng_seed = 6
rng = MersenneTwister(rng_seed)

## POMCPDPW ##

rng3 = copy(rng)
action_gen = GPSFirst(rng3)
rollout_policy = Downhill(pomdp)
tree_queries = 100
updater = ObsAdaptiveParticleFilter(pomdp, LowVarianceResampler(10_000), 0.05, rng)
node_updater = ObsAdaptiveSRFilter(pomdp, LowVarianceResampler(100), 0.05, rng)
solver = POMCPDPWSolver(next_action=action_gen,
                        tree_queries=tree_queries,
                        c = exp(35.0),
                        max_depth=mdp(pomdp).duration+2,
                        k_action=8.0,
                        alpha_action=1/20,
                        k_observation=8.0,
                        alpha_observation=1/10,
                        estimate_value=FORollout(rollout_policy),
                        # node_sr_belief_updater=node_updater,
                        rng=rng3
                       )
policy = solve(solver, pomdp)

#=
ib = initial_state_distribution(pomdp)
a = action(policy, ib)
POMCP.blink(get(policy._tree_ref))
=#

#=
using ProfileView
ib = initial_state_distribution(pomdp)
a = action(policy, ib)
Profile.clear()
@profile a = action(policy, ib)
ProfileView.view()
=#

hr = HistoryRecorder(max_steps=mdp(pomdp).duration, rng=MersenneTwister(rng_seed), show_progress=true)
hist = simulate(hr, pomdp, policy, updater)
println(discounted_reward(hist))

# using JLD
# filename = tempname()*".jld"
# JLD.save(filename, "hist", hist, "solver", solver, "pomdp", pomdp)
# println("saved to "*filename)

gr()
frames = Frames(MIME("image/png"), fps=2)
@showprogress "Creating gif..." for i in 1:n_steps(hist)
    plot(pomdp)
    push!(frames, plot!(view(hist, 1:i)))
end

filename = string(tempname(), "_pomcpdpw_powseeker.gif")
write(filename, frames)
println(filename)
run(`setsid gifview $filename`)


text:
module FNC

using LinearAlgebra, SparseArrays
using Polynomials
using NLsolve, DifferentialEquations

include("functions/chapter01.jl")
include("functions/chapter02.jl")
include("functions/chapter03.jl")
include("functions/chapter04.jl")
include("functions/chapter05.jl")
include("functions/chapter06.jl")
include("functions/chapter08.jl")
include("functions/chapter09.jl")
include("functions/chapter10.jl")
include("functions/chapter11.jl")
include("functions/chapter13.jl")

end

max_stars_repo_name: snowdj/fnc-extras

max_issues_repo_name: snowdj/fnc-extras

max_forks_repo_name: snowdj/fnc-extras


text:
# note: "contains" was manually removed from the input as it hinders parsing

struct Food
    ingredients :: Array{String, 1}
    allergens :: Array{String, 1}
end

# prints solution of both parts
function x()
    f = open("21.txt")
    N = 37
    foods = []
    for i=1:N
        l = readline(f)
        (ingr, aller) = split(l, " (")
        ing = split(ingr, " ")
        all = split(aller, ", ")
        all[end] = all[end][1:end-1]
        push!(foods, Food(ing, all))
    end

    all = Set{String}()
    for i=1:N
        for a in foods[i].allergens
            push!(all, a)
        end
    end

    sus = Set{String}()
    data = Dict{String, Set{String}}()
    for a in all
        s = Set{String}()
        f = false
        for i = 1:N
            if a in foods[i].allergens
                if f
                    t = Set{String}()
                    for x in foods[i].ingredients
                        push!(t, x)
                    end
                    intersect!(s, t)
                else
                    for x in foods[i].ingredients
                        push!(s, x)
                    end
                    f = true
                end
            end
        end
        union!(sus, s)
        data[a] = s
    end

    c = 0
    for food in foods
        for ing in food.ingredients
            if !(ing in sus)
                c += 1
            end
        end
    end
    println(c)

    arr = Array{Tuple{String,String},1}()
    i = 0
    while i < length(all)
        for (k,v) in data
            if length(v) == 1
                v0 = first(v)
                push!(arr, (k, v0))
                for (k1, v1) in data
                    delete!(v1, v0)
                end
                i += 1
            end
        end
    end
    join(map(x->x[2], sort(arr)), ",")
end


text:
mutable struct TMLEstimator <: DeterministicComposite
    Q̅::Supervised
    G::Supervised
    F::Union{LinearRegressor, LinearBinaryClassifier}
    queries::Tuple{Vararg{NamedTuple}}
    threshold::Float64
end

"""
    TMLEstimator(Q̅, G, F, query; threshold=0.005)

Implements the Targeted Minimum Loss-Based Estimator introduced by
van der Laan in https://pubmed.ncbi.nlm.nih.gov/22611591/. Two functionals of the
data generating distribution can currently be estimated:

- The classic Average Treatment Effect (ATE)
- The Interaction Average Treatment Effect (IATE) defined by Beentjes and Khamseh in
https://link.aps.org/doi/10.1103/PhysRevE.102.053314. For instance, The IATE is defined for two treatment variables as:

IATE = E[E[Y|T₁=1, T₂=1, W=w] - E[E[Y|T₁=1, T₂=0, W=w]
        - E[E[Y|T₁=0, T₂=1, W=w] + E[E[Y|T₁=0, T₂=0, W=w]

where:

- Y is the target variable (Binary)
- T = T₁, T₂ are the treatment variables (Binary)
- W are confounder variables

The TMLEstimator procedure relies on plugin estimation. Like the ATE, the IATE
requires an estimator of t,w → E[Y|T=t, W=w], an estimator of  w → p(T|w)
and an estimator of w → p(w). The empirical distribution will be used for w → p(w) all along.
The estimator of t,w → E[Y|T=t, W=w] is then fluctuated to solve the efficient influence
curve equation.

# Arguments:

- Q̅: A Supervised learner for E[Y|W, T]
- G: A Supervised learner for p(T | W)
- queries...: At least one query
- threshold: p(T | W) is truncated to this value to avoid division overflows.
"""
function TMLEstimator(Q̅::Supervised, G::Supervised, queries::Vararg{NamedTuple}; threshold=0.005::Float64)
    if Q̅ isa Probabilistic
        F = LinearBinaryClassifier(fit_intercept=false, offsetcol = :offset)
    elseif Q̅ isa Deterministic
        F = LinearRegressor(fit_intercept=false, offsetcol = :offset)
    end
    TMLEstimator(Q̅, G, F, queries, threshold)
end

###############################################################################
## Fit
###############################################################################

"""
    MLJBase.fit(tmle::TMLEstimator,
                 verbosity::Int,
                 T,
                 W,
                 y::Union{CategoricalVector{Bool}, Vector{<:Real}}

As per all MLJ inputs, T and W should respect the Tables.jl interface.
"""
function MLJBase.fit(tmle::TMLEstimator,
                 verbosity::Int,
                 T,
                 W,
                 y::Union{CategoricalVector{Bool}, Vector{<:Real}})
    Ts = source(T)
    Ws = source(W)
    ys = source(y)

    # Converting all tables to NamedTuples
    T = node(t->NamedTuple{keys(tmle.queries[1])}(Tables.columntable(t)), Ts)
    W = node(w->Tables.columntable(w), Ws)
    # intersect(keys(T), keys(W)) == [] || throw("T and W should have different column names")

    # Initial estimate of E[Y|T, W]:
    #   - The treatment variables are hot-encoded
    #   - W and T are merged
    Hmach = machine(OneHotEncoder(drop_last=true), T)
    Thot = transform(Hmach, T)

    X = node((t, w) -> merge(t, w), Thot, W)
    Q̅mach = machine(tmle.Q̅, X, ys)

    # Initial estimate of P(T|W)
    #   - T is converted to an Array
    #   - The machine is implicitely fit
    Gmach = machine(tmle.G, W, adapt(T))

    offset = compute_offset(Q̅mach, X)
    # Loop over queries that will define
    # new covariate values
    reported = []
    predicted = []
    extreme_propensity = nothing
    for (i, query) in enumerate(tmle.queries)
        indicators = indicator_fns(query)
        covariate = compute_covariate(Gmach, W, T, indicators;
                                      threshold=tmle.threshold)
        # Log extreme values
        extreme_propensity = log_over_threshold(covariate, tmle.threshold)

        # Fluctuate E[Y|T, W]
        # on the covariate and the offset
        Xfluct = fluctuation_input(covariate, offset)
        Fmach = machine(tmle.F, Xfluct, ys)

        observed_fluct = predict_mean(Fmach, Xfluct)

        queryreport = estimation_report(Fmach,
                        Q̅mach,
                        Gmach,
                        Hmach,
                        W,
                        T,
                        observed_fluct,
                        ys,
                        covariate,
                        indicators,
                        tmle.threshold,
                        query)

        push!(reported, NamedTuple{Tuple([Symbol("queryreport_$i")])}([queryreport]))
        # This is actually empty but required
        push!(predicted, observed_fluct)
    end

    predicted = hcat(predicted...)

    mach = machine(Deterministic(), Ts, Ws, ys;
            predict=predicted,
            report=(extreme_propensity_idx=extreme_propensity, merge(reported...)...)
    )
    return!(mach, tmle, verbosity)
end

max_stars_repo_name: olivierlabayle/TMLE

max_issues_repo_name: olivierlabayle/TMLE

max_forks_repo_name: olivierlabayle/TMLE


text:
# Special Pythagorean triplet

TRIPLET_SUM = 1000

function solve()
    for b in floor(Int, TRIPLET_SUM / (2 + 2^0.5)):(TRIPLET_SUM // 2)
        for a in 1:(b - 1)
            c = (a^2 + b^2)^0.5
            if a + b + c == TRIPLET_SUM
                return a * b * c
            end
        end
    end
end

print(Int(solve()))

max_stars_repo_name: wephy/project-euler

max_issues_repo_name: wephy/project-euler

max_forks_repo_name: wephy/project-euler


text:
# a simple LP example
#
#   maximize x + y
#
#   s.t. 50 x + 24 y <= 2400
#        30 x + 33 y <= 2100
#        x >= 45, y >= 5
#
#   solution: x = 45, y = 6.25, objv = 51.25

using CPLEX, Base.Test

@testset "LP 01" begin
    env = CPLEX.Env()

    # method = getparam(env, "Method")
    # println("method = $method")

    model = CPLEX.Model(env, "lp_01")
    CPLEX.set_sense!(model, :Max)

    # add variables
    CPLEX.add_var!(model, 1.0, 45., Inf)  # x
    CPLEX.add_var!(model, 1.0,  5., Inf)  # y

    # add constraints
    CPLEX.add_constr!(model, [50., 24.], '<', 2400.)
    CPLEX.add_constr!(model, [30., 33.], '<', 2100.)

    # perform optimization
    CPLEX.optimize!(model)

    sol = CPLEX.get_solution(model)
    @test sol[1] ≈ 45
    @test sol[2] ≈ 6.25

    @test CPLEX.get_objval(model) ≈ 51.25

    gc()  # test finalizers
end


text:
module C0_0
function f(q::Vector{T}, q_dot::Vector{T}, xi::T) where T
    l1 = q[1]
    l2 = q[2]
    l3 = q[3]
    l1_dot = q_dot[1]
    l2_dot = q_dot[2]
    l3_dot = q_dot[3]

    l1_dot.*(9.00239411475607e-11*l1 - 5.51654750651859e-11*l2 - 3.53674108537383e-11*l3 + 7.74272228660525e-14) - l2_dot.*(5.51654750651859e-11*l1 - 7.00120879963642e-11*l2 + 1.44764712792777e-11*l3 + 5.68275037998965e-14) - l3_dot.*(3.53674108537383e-11*l1 + 1.44764712792777e-11*l2 - 4.99826877401903e-11*l3 + 2.14973163035824e-14)
end
end


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

using Markdown
using InteractiveUtils

# ╔═╡ 22f805e8-7a7b-11eb-05e8-3dd4205997e1
begin
	using Test
	using LinearAlgebra
end

# ╔═╡ edb0b6e8-7a7a-11eb-0a88-990e1165fc1c
md"""
## \#2
"""

# ╔═╡ 357c1f92-7a7b-11eb-22bf-cd420daded2c
begin
	A = [
		1 -1 1 -1
		0  1 2  3
		0  0 2  6
		0  0 0  6
	]
	inv(A) * [
		 0
		-1//2
		-1
		-1
	]
end

# ╔═╡ d1f578ee-7a82-11eb-1375-c314d8672dd4
md"""
## \#3
"""

# ╔═╡ d5880030-7a82-11eb-3b88-a90d14844ed0
begin
	A₂ = [
		 1  1  1 1
		-1  1 -1 1
		 3  2  1 0
		 3 -2  1 0
	]
	inv(A₂) * [
		-1//6
		 0
		-1//2
		 0
	]
end

# ╔═╡ Cell order:
# ╠═22f805e8-7a7b-11eb-05e8-3dd4205997e1
# ╠═edb0b6e8-7a7a-11eb-0a88-990e1165fc1c
# ╠═357c1f92-7a7b-11eb-22bf-cd420daded2c
# ╠═d1f578ee-7a82-11eb-1375-c314d8672dd4
# ╠═d5880030-7a82-11eb-3b88-a90d14844ed0


max_stars_repo_name: AnAverageHuman/competitive

max_issues_repo_name: AnAverageHuman/competitive

max_forks_repo_name: AnAverageHuman/competitive


text:
function accuracy_thresh(ŷs, ys, thresh = 0.5)
    mean((sigmoid.(ŷs) .> thresh) .== ys)
end


text:
module CodeTools

using MacroTools

using LNR, Lazy, Requires, Compat

typealias AString AbstractString

include("utils.jl")
include("eval.jl")
include("module.jl")
include("summaries.jl")
include("completions.jl")
include("doc.jl")

end # module


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

mutable struct IoK8sApiAutoscalingV2beta2ObjectMetricSource <: SwaggerModel
    describedObject::Any # spec type: Union{ Nothing, IoK8sApiAutoscalingV2beta2CrossVersionObjectReference } # spec name: describedObject
    metric::Any # spec type: Union{ Nothing, IoK8sApiAutoscalingV2beta2MetricIdentifier } # spec name: metric
    target::Any # spec type: Union{ Nothing, IoK8sApiAutoscalingV2beta2MetricTarget } # spec name: target

    function IoK8sApiAutoscalingV2beta2ObjectMetricSource(;describedObject=nothing, metric=nothing, target=nothing)
        o = new()
        validate_property(IoK8sApiAutoscalingV2beta2ObjectMetricSource, Symbol("describedObject"), describedObject)
        setfield!(o, Symbol("describedObject"), describedObject)
        validate_property(IoK8sApiAutoscalingV2beta2ObjectMetricSource, Symbol("metric"), metric)
        setfield!(o, Symbol("metric"), metric)
        validate_property(IoK8sApiAutoscalingV2beta2ObjectMetricSource, Symbol("target"), target)
        setfield!(o, Symbol("target"), target)
        o
    end
end # type IoK8sApiAutoscalingV2beta2ObjectMetricSource

const _property_map_IoK8sApiAutoscalingV2beta2ObjectMetricSource = Dict{Symbol,Symbol}(Symbol("describedObject")=>Symbol("describedObject"), Symbol("metric")=>Symbol("metric"), Symbol("target")=>Symbol("target"))
const _property_types_IoK8sApiAutoscalingV2beta2ObjectMetricSource = Dict{Symbol,String}(Symbol("describedObject")=>"IoK8sApiAutoscalingV2beta2CrossVersionObjectReference", Symbol("metric")=>"IoK8sApiAutoscalingV2beta2MetricIdentifier", Symbol("target")=>"IoK8sApiAutoscalingV2beta2MetricTarget")
Base.propertynames(::Type{ IoK8sApiAutoscalingV2beta2ObjectMetricSource }) = collect(keys(_property_map_IoK8sApiAutoscalingV2beta2ObjectMetricSource))
Swagger.property_type(::Type{ IoK8sApiAutoscalingV2beta2ObjectMetricSource }, name::Symbol) = Union{Nothing,eval(Meta.parse(_property_types_IoK8sApiAutoscalingV2beta2ObjectMetricSource[name]))}
Swagger.field_name(::Type{ IoK8sApiAutoscalingV2beta2ObjectMetricSource }, property_name::Symbol) =  _property_map_IoK8sApiAutoscalingV2beta2ObjectMetricSource[property_name]

function check_required(o::IoK8sApiAutoscalingV2beta2ObjectMetricSource)
    (getproperty(o, Symbol("describedObject")) === nothing) && (return false)
    (getproperty(o, Symbol("metric")) === nothing) && (return false)
    (getproperty(o, Symbol("target")) === nothing) && (return false)
    true
end

function validate_property(::Type{ IoK8sApiAutoscalingV2beta2ObjectMetricSource }, name::Symbol, val)
end