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

using Markdown
using InteractiveUtils

# This Pluto notebook uses @bind for interactivity. When running this notebook outside of Pluto, the following 'mock version' of @bind gives bound variables a default value (instead of an error).
macro bind(def, element)
    quote
        local el = $(esc(element))
        global $(esc(def)) = Core.applicable(Base.get, el) ? Base.get(el) : missing
        el
    end
end

# ╔═╡ 10aba104-a93f-11eb-29ed-4136efb69580
begin
    import Pkg
    Pkg.activate(mktempdir())
    Pkg.add([
        Pkg.PackageSpec(name="Latexify"),
        Pkg.PackageSpec(name="LaTeXStrings"),
        Pkg.PackageSpec(name="Plots"),
		Pkg.PackageSpec(name="ColorSchemes"),
		Pkg.PackageSpec(name="PlotThemes"),
		Pkg.PackageSpec(name="SymEngine"),
		Pkg.PackageSpec(name="PlutoUI")
    ])
    using Latexify
    using LaTeXStrings
	using Plots
	using ColorSchemes
	using PlotThemes
	using SymEngine
	using PlutoUI
end

# ╔═╡ 1a671d44-ac1c-4871-8bcd-7d16642f34e0
begin
	Pkg.add(url="https://github.com/Pocket-titan/DarkMode")
    import DarkMode
end

# ╔═╡ f5b9db94-6ab4-4abb-ac29-6b03ac4dcb16
begin
	md"# Filtro Variables de Estado

Tomamos el siguiente circuito, con el objetivo de analizar y simular sus transferencias a V1, V2, V3."
end

# ╔═╡ 1c78f0c7-8485-416f-8789-716571d875b6
md"Planteando las ecuaciones de nodos, se llega al siguiente sistema de 4 ecuaciones con 4 incógnitas y Vi parámetro"

# ╔═╡ 9978620e-7bcc-4c13-b3d9-0dc714e95dae
begin
	nodo1 = :(1/R_4*(V_1-V_X)+1/R_3*(V_3-V_X) = 0)
	nodo2 = :(1/R_1*V_1 + jwC_2*V_3 = 0)
	nodo3 = :(1/R_2*V_2 + jwC_3*V_3 = 0)
	nodo4 = :(1/R_5*(V_i-V_X) + 1/R_6*(V_2-V_X) = 0)
	latexraw(nodo2)
	L"$\begin{aligned}
	\tfrac{1}{R_{4}} \cdot \left( V_{1} - V_{X} \right) + \tfrac{1}{R_{3}} \cdot \left( V_{3} - V_{X} \right) &= 0\\
	\tfrac{1}{R_{1}} \cdot V_{1} + jwC_{2} \cdot V_{2} &= 0\\
	\tfrac{1}{R_{2}} \cdot V_{2} + jwC_{3} \cdot V_{3} &= 0\\
	\tfrac{1}{R_{5}} \cdot \left( V_{i} - V_{X} \right) + \tfrac{1}{R_{6}} \cdot \left( V_{2} - V_{X} \right) &= 0
	\end{aligned}$"
end

# ╔═╡ 7c4cd737-ea86-4d02-92bf-774f8fb04ad0
md"Expresado matricialmente:"

# ╔═╡ 248a48d7-493b-4506-ae69-3fdbbf82f7c9
L"$\begin{bmatrix}
	\tfrac{1}{R_{4}} + \tfrac{1}{R_3} & -\tfrac{1}{R_4} & 0 & -\tfrac{1}{R_3}\\
	0 & \tfrac{1}{R_1} & jwC_2 & 0\\
	0 & 0 & \tfrac{1}{R_2} & jwC_3\\
	\tfrac{1}{R_{5}} + \tfrac{1}{R_5} & 0 & -\tfrac{1}{R_6} & 0\\
\end{bmatrix}
\begin{bmatrix}
V_X\\ V_2\\ V_3\\ V_4\\
\end{bmatrix}
=
\begin{bmatrix}
0\\ 0\\ 0\\ \tfrac{1}{R_5}\\
\end{bmatrix} V_i$"

# ╔═╡ 44009b12-cbe2-4d18-ba4a-9b4bcf71f465
s = symbols(:s)

# ╔═╡ c5d3db4a-7931-49ae-86dc-a8db73ae2ba1
md"(La funcionalidad interactiva de Pluto Notebooks es todavía incompatible con Github Pages)"

# ╔═╡ fe9697e4-135a-4dd2-9bf1-3d47f6a9da0e
#Rs = @bind R Slider(1:100; default=1, show_value=true)
Rs = 1

# ╔═╡ e01e76cf-6384-4534-8382-eea9a040e25d
#Cs = @bind C Slider(1:100; default=100, show_value=true)
Cs = 100

# ╔═╡ f25d8bcd-5887-4332-b54e-d9fb540825f2
begin
	R1 = Rs*1e3;
	R2 = Rs*1e3;
	R3 = Rs*1e3;
	R4 = Rs*1e3;
	R5 = Rs*1e3;
	R6 = Rs*1e3;
	C2 = Cs*1e-9;
	C3 = Cs*1e-9;
	md"Usando el paquete SymEngine.jl, se puede llegar a una resolución numérica-simbólica, definiendo la variable simbólica s y los valores numéricos de R y C. Los selectores permitirían asignar valores a R en kΩ y C en nF para observar su efecto en la transferencia."
end

# ╔═╡ 71d5083b-ba6c-4a77-a0c9-3b79a282e636
md"En función de eso, se instancia la matriz de admitancias Y y el vector de parámetros i."

# ╔═╡ 1092d595-c1db-4a08-9d18-dacc808e8ae9
Y = [1/R4+1/R3 -1/R4 0    -1/R3;
     0          1/R1 s*C2  0;
     0          0    1/R2  s*C3;
     1/R5+1/R5  0   -1/R6  0]

# ╔═╡ e73c6a27-f6ec-4e05-8fc2-5e277a6c05ce
i = [0; 0; 0; 1/R3]

# ╔═╡ 5a7e24aa-fefb-4e04-8440-2598d3004165
begin
	H = expand.(inv(Y)*i)
	md"SymEngine.jl permite resolver el sistema YH = i semi-simbólicamente, obteniendo así el vector de transferencias H = [Hx; H1; H2; H3] = inv(Y) i, de donde se pueden extraer las transferencias."
end

# ╔═╡ 2324a9f8-a187-4435-81e8-df2c4b44eb90
begin
	H1(s) = H[2]
	H2(s) = H[3]
	H3(s) = H[4]
end

# ╔═╡ 4b87aa13-e337-4512-b477-8e125da89e39
md"Las transferencias corresponden respectivamente a H1 pasa altos, H2 pasa banda, y H3 pasa bajos."

# ╔═╡ eb0ea7cc-49f9-4e15-9382-193b52ad485c
latexstring("H_1(s) ="*latexraw(H1(s)))

# ╔═╡ 8299fea2-4f51-45b9-a98a-1fb11b511d80
latexstring("H_2(s) ="*latexraw(H2(s)))

# ╔═╡ 1f845ecc-353a-42c8-80dc-159bed6d1652
latexstring("H_3(s) ="*latexraw(H3(s)))

# ╔═╡ 5e678738-638b-44c3-9d6a-df86c3788e01
md"Definiendo la variable simbólica w, SymEngine.jl permite sustituir s = jw en la transferencia H(s) para visualizar las transferencias como diagrama de Bode."

# ╔═╡ eea7c707-8c2f-48e4-86d9-4ede8b111fab
begin
	w = symbols(:w)
	H1w(w) = subs(H1(s), s => w*im)
	H2w(w) = subs(H2(s), s => w*im)
	H3w(w) = subs(H3(s), s => w*im)
end

# ╔═╡ b0fec057-d05b-4a9b-835e-cd267465ff0a
md"A continuación se presentan las gráficas de |H(w)|. Se verifica que H1 corresponde a la salida pasa altos, H2 a la salida pasa banda, y H3 a la salida pasa bajos."

# ╔═╡ 44bac19d-b875-478e-8a68-2b640765f8cf
begin
	theme(:juno)
	pH1 = plot(abs(H1w(w)), 1, 1e8, xaxis=:log, yaxis=:log, label="H1")
	pH2 = plot(abs(H2w(w)), 1, 1e8, xaxis=:log, yaxis=:log, label="H2")
	pH3 = plot(abs(H3w(w)), 1, 1e8, xaxis=:log, yaxis=:log, label="H3", xlabel="w")
	plot(pH1, pH2, pH3, layout = (3,1))
end

# ╔═╡ 1df8580c-d92d-4c17-b9f7-d934379c80ad
md"Ahora se procede a una resolución puramente simbólica de la transferencia del filtro, definiendo R y C como variables. Por simplicidad del resultado final suponemos resistencias y capacitores de igual valor."

# ╔═╡ cfee3c91-9ec7-4e1f-9fda-f7b9fb2407e0
@vars R C

# ╔═╡ 51bee854-294a-44c3-89ae-392969e53c08
Ys = [1/R+1/R -1/R  0   -1/R;
      0        1/R  s*C  0;
      0        0    1/R  s*C;
     1 /R+1/R  0   -1/R  0]

# ╔═╡ 690f4c2c-a454-4c18-b798-13f4e2c08eb8
is = [0; 0; 0; 1/R]

# ╔═╡ 8d914ac6-b6fd-4ab9-a7d4-40645205c3ea
begin
	Hs = expand.(inv(Ys)*is)
	Hs1(s) = Hs[2]
	Hs2(s) = Hs[3]
	Hs3(s) = Hs[4]
	md"Con un procedimiento idéntico al anterior se puede obtener una expresión puramente simbólica de las transferencias H1, H2, H3:"
end

# ╔═╡ 0d8fa251-0e1e-4f24-9e4b-5eab54cb483d
latexstring("H_1(s) ="*latexraw(Hs1(s)))

# ╔═╡ 3ecde8b1-5b0c-4a78-8768-46a5cdf25dff
latexstring("H_2(s) ="*latexraw(Hs2(s)))

# ╔═╡ 28fe2ae9-9f71-4dfb-a8a0-bcb925498978
latexstring("H_3(s) ="*latexraw(Hs3(s)))

# ╔═╡ 35ee3ea7-465e-4890-8bd9-98066ddfb10a
md"O bien, escrito de la forma habitual:"

# ╔═╡ 7642f158-b9ab-4ca5-af98-0ca8fb9c053c
	L"$\begin{aligned}
	H_1(s) &= \frac{R^2C^2s^2}{R^2C^2s+RCs+1}\\[1em]
	H_2(s) &= \frac{RCs}{R^2C^2s+RCs+1}\\[1em]
	H_1(s) &= \frac{1}{R^2C^2s+RCs+1}\\
	\end{aligned}$"

# ╔═╡ 8a097b97-65bb-4b9f-87c2-25226ab2a7fc
md"Es decir, resultan tres filtros con frecuencia de corte wc=RC, que en el caso particular correspondiente a la resolución numérica de R = 1kΩ y C = 100nF retorna una frecuencia wc = 10000 rad/s, consistente con lo visto en los diagramas de Bode."

# ╔═╡ 6ec8a49a-fb90-493a-8301-9aa23c482c25
begin
	struct Wow
		filename
	end

	function Base.show(io::IO, ::MIME"image/png", w::Wow)
		write(io, read(w.filename))
	end
end

# ╔═╡ 98c4b949-7464-4690-92cc-93e205a35072
Wow("Images/circuito.png")

# ╔═╡ e8fca5e5-2263-4e73-9481-292a528961c1
DarkMode.enable(theme="nord")

# ╔═╡ Cell order:
# ╟─f5b9db94-6ab4-4abb-ac29-6b03ac4dcb16
# ╟─98c4b949-7464-4690-92cc-93e205a35072
# ╟─1c78f0c7-8485-416f-8789-716571d875b6
# ╟─9978620e-7bcc-4c13-b3d9-0dc714e95dae
# ╟─7c4cd737-ea86-4d02-92bf-774f8fb04ad0
# ╟─248a48d7-493b-4506-ae69-3fdbbf82f7c9
# ╟─f25d8bcd-5887-4332-b54e-d9fb540825f2
# ╠═44009b12-cbe2-4d18-ba4a-9b4bcf71f465
# ╟─c5d3db4a-7931-49ae-86dc-a8db73ae2ba1
# ╠═fe9697e4-135a-4dd2-9bf1-3d47f6a9da0e
# ╠═e01e76cf-6384-4534-8382-eea9a040e25d
# ╟─71d5083b-ba6c-4a77-a0c9-3b79a282e636
# ╠═1092d595-c1db-4a08-9d18-dacc808e8ae9
# ╠═e73c6a27-f6ec-4e05-8fc2-5e277a6c05ce
# ╟─5a7e24aa-fefb-4e04-8440-2598d3004165
# ╠═2324a9f8-a187-4435-81e8-df2c4b44eb90
# ╟─4b87aa13-e337-4512-b477-8e125da89e39
# ╠═eb0ea7cc-49f9-4e15-9382-193b52ad485c
# ╠═8299fea2-4f51-45b9-a98a-1fb11b511d80
# ╠═1f845ecc-353a-42c8-80dc-159bed6d1652
# ╟─5e678738-638b-44c3-9d6a-df86c3788e01
# ╠═eea7c707-8c2f-48e4-86d9-4ede8b111fab
# ╟─b0fec057-d05b-4a9b-835e-cd267465ff0a
# ╟─44bac19d-b875-478e-8a68-2b640765f8cf
# ╟─1df8580c-d92d-4c17-b9f7-d934379c80ad
# ╟─cfee3c91-9ec7-4e1f-9fda-f7b9fb2407e0
# ╟─51bee854-294a-44c3-89ae-392969e53c08
# ╟─690f4c2c-a454-4c18-b798-13f4e2c08eb8
# ╟─8d914ac6-b6fd-4ab9-a7d4-40645205c3ea
# ╟─0d8fa251-0e1e-4f24-9e4b-5eab54cb483d
# ╟─3ecde8b1-5b0c-4a78-8768-46a5cdf25dff
# ╟─28fe2ae9-9f71-4dfb-a8a0-bcb925498978
# ╟─35ee3ea7-465e-4890-8bd9-98066ddfb10a
# ╟─7642f158-b9ab-4ca5-af98-0ca8fb9c053c
# ╟─8a097b97-65bb-4b9f-87c2-25226ab2a7fc
# ╟─6ec8a49a-fb90-493a-8301-9aa23c482c25
# ╟─10aba104-a93f-11eb-29ed-4136efb69580
# ╟─1a671d44-ac1c-4871-8bcd-7d16642f34e0
# ╟─e8fca5e5-2263-4e73-9481-292a528961c1

max_stars_repo_path: 1 - Filtro Variables Estado/notebook.jl

max_stars_repo_name: Matrsz/Teoria-Circuitos-II-UTN

max_issues_repo_path: 1 - Filtro Variables Estado/notebook.jl

max_issues_repo_name: Matrsz/Teoria-Circuitos-II-UTN

max_forks_repo_path: 1 - Filtro Variables Estado/notebook.jl

max_forks_repo_name: Matrsz/Teoria-Circuitos-II-UTN


text:
# Timers.

"""
Global vector of timers `(deadline, condition)` sorted by `deadline`.
"""
struct Timer
    deadline::Float64
    f::Any
end

Base.isless(a::Timer, b::Timer) = isless(a.deadline, b.deadline)

const timer_vector = Timer[]
const timer_lock = Threads.SpinLock()

function dump_timer_state(msg="")
    now = time()
    printerr("Timers $msg: ", [round(t.deadline - now; digits=2)
                               for t in timer_vector])
end

"""
    UnixIO.sleep(seconds)

Block the current task for a specified number of seconds.
(Waits for a timer based on `UnixIO.poll_wait()`).
"""
sleep(seconds::Float64) = wait_timeout(sleep_condition, seconds)

const sleep_condition = Base.ThreadSynchronizer()

"""
Restart the clock source responsible for calling `notify_timers()`.
The clock source will be sleeping based on `next_timer_deadline_ms()`,
if timer is registered with a nearer deadline, the clock source needs
to be restarted (it will call `next_timer_deadline_ms()` again and go
back to sleep).
"""
@db 3 function restart_timer_clock()
    wakeup_poll(poll_queue)
end

"""
    next_timer_deadline_ms(timeout_ms)

Time from now to the next timer deadline in ms.
If no timers due before `timeout_ms` return `timeout_ms`.
"""
function next_timer_deadline_ms(timeout_ms::Int)
    @dblock timer_lock begin
        deadline = next_timer_deadline()
        dt = deadline - time()
        if dt < (timeout_ms/1000)
            timeout_ms = max(0, round(Int, dt * 1000))
        end
    end
    return timeout_ms
end

function next_timer_deadline()
    assert_havelock(timer_lock)
    isempty(timer_vector) ? Inf : timer_vector[1].deadline
end

const inf_timer = Timer(Inf, nothing)

"""
    register_timer(f, deadline)

Register `f`` to be called at `time()` `deadline`.
"""
@db 5 function register_timer(f, deadline)
    if deadline == Inf
        return inf_timer
    end
    @db 5 return register_timer(Timer(deadline, f))
end

@db 3 function register_timer(t::Timer)
    i = 0
    @dblock timer_lock begin
        i = searchsortedfirst(timer_vector, t)
        insert!(timer_vector, i, t)
    end
    if i == 1                                    ;@db 3 "restart_timer_clock()"
        restart_timer_clock()
    end
    @db 3 return t
end

@db 3 function cancel_timer(t::Timer)
    @require t.deadline !=  Inf
    @dblock timer_lock begin
        if t.deadline >= next_timer_deadline()
            i = searchsortedfirst(timer_vector, t)
            if timer_vector[i] == t                              ;@db 3 "⏱ 🚫"
                deleteat!(timer_vector, i)
            end
        end
    end
end

"""
Cancel a timer `(deadline, condition)` without notifying `condition`.
"""
@db 5 function Base.close(t::Timer)
    if t.deadline != Inf
        cancel_timer(t)
    end
    nothing
end

"""
Notify timers whos `deadline` has passed and remove from vector.
"""
@db 6 function notify_timers(now = time())
    ready=[]
    @dblock timer_lock begin
        expired = 0
        for (i, t) in enumerate(timer_vector)
            if t.deadline > now
                break
            end
            push!(ready, t.f)
            expired = i
        end
        deleteat!(timer_vector, 1:expired)
    end
    for f in ready                                            ;@db 3 "⏱ -> $f"
        f()
    end
    nothing
end

"""
    wait_until(condition, [predicate=()->false], deadline)

If `predicate()` is false, repeatently `wait(condition)` until
`predicate()` is true (or until `deadline` is reached).
"""
@db 3 function wait_until(condition::Base.GenericCondition,
                          predicate,
                          deadline::Float64)
    if predicate()
        return
    end
    @dblock condition begin
        timer = register_timer(deadline) do
            @dblock condition notify(condition)
        end
        try
            while !predicate() && time() < deadline         ;@db 3 "waiting..."
                wait(condition)
            end
        finally
            close(timer)
        end
    end
    nothing
end

wait_until(c, d) = wait_until(c, ()->false, d)
@selfdoc wait_timeout(c, t) = wait_until(c, time() + Float64(t))

# Pretty Printing.

function dbshow(io::IO, t::Timer)
    print(io, "Timer(", debug_time(t.deadline),
                        " [", round(t.deadline - time(); digits=2), "])")
end

# End of file: timer.jl

ext: jl


text:
function GridLayout(spec::GridLayoutSpec)
    gl = GridLayout(; spec.kwargs...)

    for ((rows, cols, side), content) in spec.content
        if content isa GridLayoutSpec
            content = GridLayout(content)
        end
        gl[rows, cols, side] = content
    end

    gl
end

to_valid_pos(pos::Tuple{Any, Any}) = (pos[1], pos[2], Inner())
to_valid_pos(pos::Tuple{Any, Any, Any}) = pos

function GridLayoutSpec(content::Vector{<:Pair}; kwargs...)
    spec_content::Vector{Pair{Tuple{Indexables, Indexables, Side}, Any}} =
        map(content) do (pos, content)
            validpos = to_valid_pos(pos)
            return validpos => content
        end
    GridLayoutSpec(spec_content, Dict{Symbol, Any}(kwargs))
end


text:
// Assigning int to double variable.

int main () {
 double x = 1;
 return 0 ;
}

max_stars_repo_name: jwilk/jtc

max_issues_repo_name: jwilk/jtc

max_forks_repo_name: jwilk/jtc


text:
# NOTE: This code is purely to test the coverage of the JSON to test converter.
# It prints the number of exercises that the script works on.

using JSON, HTTP

# Include the code we want to test the coverage of
include("JSONtoTests.jl")

# Print the slugs of the exercises that the code fails on.
println("The converter fails on the following exercises: ")

# Initialize an error counter to keep track of the number of exercises our script does not work on.
error_count = 0

# Initialize a counter for keeping track of how many exercises don't have canonical data in problem-specifications.
not_existing = 0

# Get all exercises and iterate over them

slugs = collect(readdir("../exercises"))

# Comment out if you want to run on all exercises in problem-specifications
# slugs = ["accumulate","acronym","affine-cipher","all-your-base","allergies","alphametics","anagram","armstrong-numbers","atbash-cipher","bank-account","beer-song","binary","binary-search","binary-search-tree","bob","book-store","bowling","change","circular-buffer","clock","collatz-conjecture","complex-numbers","connect","counter","crypto-square","custom-set","darts","diamond","difference-of-squares","diffie-hellman","dnd-character","dominoes","dot-dsl","error-handling","etl","flatten-array","food-chain","forth","gigasecond","go-counting","grade-school","grains","grep","hamming","hangman","hello-world","hexadecimal","high-scores","house","isbn-verifier","isogram","kindergarten-garden","knapsack","largest-series-product","leap","ledger","lens-person","linked-list","list-ops","luhn","markdown","matching-brackets","matrix","meetup","micro-blog","minesweeper","nth-prime","nucleotide-codons","nucleotide-count","ocr-numbers","octal","paasio","palindrome-products","pangram","parallel-letter-frequency","pascals-triangle","perfect-numbers","phone-number","pig-latin","point-mutations","poker","pov","prime-factors","protein-translation","proverb","pythagorean-triplet","queen-attack","rail-fence-cipher","raindrops","rational-numbers","react","rectangles","resistor-color","resistor-color-duo","resistor-color-trio","rest-api","reverse-string","rna-transcription","robot-name","robot-simulator","roman-numerals","rotational-cipher","run-length-encoding","saddle-points","satellite","say","scale-generator","scrabble-score","secret-handshake","series","sgf-parsing","sieve","simple-cipher","simple-linked-list","space-age","spiral-matrix","strain","sublist","sum-of-multiples","tournament","transpose","tree-building","triangle","trinary","twelve-days","two-bucket","two-fer","variable-length-quantity","word-count","word-search","wordy","yacht","zebra-puzzle","zipper"]

for i in slugs
    # Make error_count global to change its value in the try-catch block.
    global error_count
    global not_existing

    # Try to run the code on the exercise.
    try
        main(i)
    catch error
        # If exercise does not exist in problem-specifications repo, increment number of exercies
        # that are not on problem-specifications and move on.
        if isa(error,NotFoundException)
            not_existing += 1
            continue
        end

        # If code doesn't work, then we can print the slug of the exercise. (slug is indented, and bulleted)
        println("-   $i")

        # If code doesn't work, increment the error count, and continue to next exercise
        error_count += 1
        continue
    end
end

# Print final number of passes and fails.
println("Pass: $(length(slugs) - error_count - not_existing), Fail = $error_count, Neutral = $not_existing")

max_stars_repo_name: soumitradev/julia

max_issues_repo_name: soumitradev/julia

max_forks_repo_name: soumitradev/julia


text:
# Author: Ali Siahkoohi, alisk@gatech.edu
# Date: February 2021
# Copyright: Georgia Institute of Technology, 2021

export Mask

mutable struct Mask
    M::AbstractArray{Float32,4}
    fac::Float32
end

function Mask(shape::Array{Int64,1}, fac::Float32)

    nzeros = Int(prod(shape[1:2]) * fac)

    M = ones(Float32, prod(shape[1:2]))
    M[sort(randperm(prod(shape[1:2]))[1:nzeros])] .= 0f0
    M = reshape(M, (shape[1], shape[2], 1, 1))

    return Mask(M, fac)
end

function *(Mask::Mask , X::AbstractArray{Float32,4})
    return Mask.M .* X
end


text:
# Model
include_model("hopper")
model_ft = free_time_model(no_slip_model(model))

function fd(model::Hopper, x⁺, x, u, w, h, t)
	q3 = view(x⁺, model.nq .+ (1:model.nq))
	q2⁺ = view(x⁺, 1:model.nq)
	q2⁻ = view(x, model.nq .+ (1:model.nq))
	q1 = view(x, 1:model.nq)
	u_ctrl = view(u, model.idx_u)
	λ = view(u, model.idx_λ)
	b = view(u, model.idx_b)
	h = u[end]

	[q2⁺ - q2⁻;
	((1.0 / h) * (M_func(model, q1) * (SVector{4}(q2⁺) - SVector{4}(q1))
	- M_func(model, q2⁺) * (SVector{4}(q3) - SVector{4}(q2⁺)))
	+ transpose(B_func(model, q3)) * SVector{2}(u_ctrl)
	+ transpose(N_func(model, q3)) * SVector{1}(λ)
	+ transpose(P_func(model, q3)) * SVector{2}(b)
	- h * G_func(model, q2⁺))]
end

function no_slip(model::Hopper, x⁺, u, h)
	q3 = view(x⁺, model.nq .+ (1:model.nq))
	q2 = view(x⁺, 1:model.nq)
	λ = view(u, model.idx_λ)
	s = view(u, model.idx_s)

	h = u[end]
	return s[1] - (λ' * _P_func(model, q3) * (q3 - q2) / h)[1]
end

# Horizon
T = 21

# Time step
tf = 1.0
h = tf / (T - 1)

# Bounds
_uu = Inf * ones(model_ft.m)
_uu[model_ft.idx_u] .= 25.0
_uu[end] = 5.0 * h
_ul = zeros(model_ft.m)
_ul[model_ft.idx_u] .= -25.0
_ul[end] = 0.1 * h
ul, uu = control_bounds(model_ft, T, _ul, _uu)

# Initial and final states
q1 = [0.0, 0.5 , 0.0, 0.5]
q_right = [-0.25, 0.5 + 0.25, pi / 2.0, 0.25]
q_top = [-0.5, 0.5 + 0.5, pi, 0.25]
q_left = [-0.75, 0.5 + 0.25, 3.0 * pi / 2.0, 0.25]
qT = [-1.0, 0.5,  2.0 * pi, 0.5]

xl, xu = state_bounds(model_ft, T,
		[model_ft.qL; model_ft.qL],
		[model_ft.qU; model_ft.qU],
        x1 = [q1; q1],
		xT = [Inf * ones(model_ft.nq); qT])

q_ref = [linear_interpolation(q1, q_right, 6)...,
         linear_interpolation(q_right, q_top, 6)[2:end]...,
         linear_interpolation(q_top, q_left, 6)[2:end]...,
         linear_interpolation(q_left, qT, 6)[2:end]...]

x_ref = configuration_to_state(q_ref)

# Objective
include_objective("velocity")
obj_tracking = quadratic_time_tracking_objective(
    [Diagonal(zeros(model_ft.n)) for t = 1:T],
    [Diagonal([1.0e-1, 1.0e-1, zeros(model_ft.m - model_ft.nu)...]) for t = 1:T-1],
    [zeros(model_ft.n) for t = 1:T],
    [zeros(model_ft.m) for t = 1:T-1],
    1.0)

obj_contact_penalty = PenaltyObjective(1.0e3, model_ft.m - 1)
obj_velocity = velocity_objective([Diagonal(ones(model_ft.nq)) for t = 1:T],
	model_ft.nq)
obj = MultiObjective([obj_tracking, obj_velocity, obj_contact_penalty])

# Constraints
include_constraints(["free_time", "contact_no_slip"])
con_free_time = free_time_constraints(T)
con_contact = contact_no_slip_constraints(model_ft, T)
con = multiple_constraints([con_free_time, con_contact])

# Problem
prob = trajectory_optimization_problem(
			   model_ft,
               obj,
               T,
               xl = xl,
               xu = xu,
               ul = ul,
               uu = uu,
               con = con)

# Trajectory initialization
x0 = deepcopy(x_ref) # linear interpolation on state
u0 = [[1.0e-3 * rand(model_ft.m-1); h] for t = 1:T-1] # random controls

# Pack trajectories into vector
z0 = pack(x0, u0, prob)

#NOTE: may need to run examples multiple times to get good trajectories
# Solve nominal problem
include_snopt()

@time z̄ , info = solve(prob, copy(z0),
	nlp = :SNOPT7,
	tol = 1.0e-3, c_tol = 1.0e-3, mapl = 5)

x̄, ū = unpack(z̄, prob)

@show ū[4][end]
@show check_slack(z̄, prob)

using Plots
tf, t, h = get_time(ū)
plot(t[1:end-1], hcat(ū...)[1:2,:]', linetype=:steppost,
	xlabel="time (s)", ylabel = "control",
	label = ["angle" "length"],
	width = 2.0, legend = :top)
plot(t[1:end-1], h, linetype=:steppost)

include(joinpath(pwd(), "models/visualize.jl"))
vis = Visualizer()
render(vis)
visualize!(vis, model_ft, state_to_configuration(x̄), Δt = h[1])


text:
using Test
using ParametricOptInterface
using MathOptInterface
using GLPK
using Ipopt
using JuMP

const POI = ParametricOptInterface
const MOI = MathOptInterface
const MOIU = MathOptInterface.Utilities

const ATOL = 1e-4

include("production_problem_test.jl")
include("basic_tests.jl")
include("quad_tests.jl")
include("jump_tests.jl")


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

using Markdown
using InteractiveUtils

# This Pluto notebook uses @bind for interactivity. When running this notebook outside of Pluto, the following 'mock version' of @bind gives bound variables a default value (instead of an error).
macro bind(def, element)
    quote
        local el = $(esc(element))
        global $(esc(def)) = Core.applicable(Base.get, el) ? Base.get(el) : missing
        el
    end
end

# ╔═╡ f9509a30-3a3c-11eb-16b2-6d2d967af99f
begin
	using Pkg
	Pkg.add(["FFTW","Plots","CSV","DataFrames","Statistics","PlutoUI","Peaks"])
	using FFTW, Plots, CSV, DataFrames, Statistics, PlutoUI, Peaks
end

# ╔═╡ 39d10dae-3a3d-11eb-13ea-8125fd517e4e
data = CSV.read("data/stream_simple.csv",DataFrame)

# ╔═╡ 91202802-3afe-11eb-0bec-27e17fc805fc
plotly()

# ╔═╡ ea0ad14e-3af8-11eb-2b1c-91c08a99d77d
function runave(data;alpha=0.9)
	ave = mean(data)
	runs = []
	for i in 1:length(data)
		ave = ave * alpha + (1-alpha) * data[i]
		push!(runs,ave)
	end
	return runs
end

# ╔═╡ 9f688650-3e2d-11eb-0d8d-c1fed0e55490
function peaktime(data)
	res = []
	for i in 2:length(data)
		if data[i] - data[i-1] > 100
			push!(res,50)
		else
			push!(res,data[i] - data[i-1])
		end
	end
	return res
end

# ╔═╡ 9c9c1460-3b02-11eb-0631-3d37cd287381
@bind α Slider(0.8:0.01:0.99;default=0.9)

# ╔═╡ 52643a20-3af9-11eb-13ea-91dc64e1ff1e
run_ave = Float64.(runave(data.Count;alpha=α))

# ╔═╡ f9da7cd0-3e30-11eb-2c26-6574a5f1ec7f
p = plot(data.Time,run_ave);

# ╔═╡ aba20fa0-3e2c-11eb-1115-61101d91e01f
peaks = peakprom(run_ave,20)[1]

# ╔═╡ 39381b30-3e31-11eb-0950-a30f253095af
scatter!(p,data.Time[peaks],run_ave[peaks])

# ╔═╡ d22a4430-3e3b-11eb-1e0b-951e5b602f99
peaks[peakprom(run_ave,20)[2] .> 0.2]

# ╔═╡ 8bf03660-3e30-11eb-188b-e9e1bdd19a73
res = peaktime(peaks)

# ╔═╡ fe864a40-3e33-11eb-1b29-5905138e815f
bpm = 60.0 ./ (res ./ 10)

# ╔═╡ fabb2030-3e2d-11eb-00b5-2b878c128f8b
plot(data.Time[peaks[2:end]],[bpm,runave(bpm)])

# ╔═╡ Cell order:
# ╠═f9509a30-3a3c-11eb-16b2-6d2d967af99f
# ╠═39d10dae-3a3d-11eb-13ea-8125fd517e4e
# ╠═91202802-3afe-11eb-0bec-27e17fc805fc
# ╠═ea0ad14e-3af8-11eb-2b1c-91c08a99d77d
# ╠═9f688650-3e2d-11eb-0d8d-c1fed0e55490
# ╠═52643a20-3af9-11eb-13ea-91dc64e1ff1e
# ╠═9c9c1460-3b02-11eb-0631-3d37cd287381
# ╠═f9da7cd0-3e30-11eb-2c26-6574a5f1ec7f
# ╠═39381b30-3e31-11eb-0950-a30f253095af
# ╠═aba20fa0-3e2c-11eb-1115-61101d91e01f
# ╠═d22a4430-3e3b-11eb-1e0b-951e5b602f99
# ╠═8bf03660-3e30-11eb-188b-e9e1bdd19a73
# ╠═fe864a40-3e33-11eb-1b29-5905138e815f
# ╠═fabb2030-3e2d-11eb-00b5-2b878c128f8b

max_stars_repo_name: RobViren/SleepWyzely

max_issues_repo_name: RobViren/SleepWyzely

max_forks_repo_name: RobViren/SleepWyzely


text:
function compute_obs(O::Thermo, S::Array{Int, 3}, M::Model)
    L = M.dims[1]
    P = length(O.D.∂E[:α])
    m = P + 1
    E, Z, X, ZZ, mag = 0.0, 0.0, 0.0, 0.0, 0.0
    for var in [:α, :β, :γ]
        O.D.tmp_∂E[var] .= zeros(P)
    end
    O.D.tmp_X_∂X .= zeros(P)

    for j in 1:L, i in 1:L
        if i < L && j < L
            zz = S[i, j, m] * S[i+1, j, m]
            zz += S[i, j, m] * S[i, j+1, m]
        elseif j == L && i < L
            zz = S[i, j, m] * S[i+1, j, m]
        elseif i == L && j < L
            zz = S[i, j, m] * S[i, j+1, m]
        else
            zz = 0
        end
        x = exp(-2.0 * M.J[:X][m] * S[i, j, m] * S[i, j, m+1]) # SLOW LINE 3
        # z = S[i, j, m]
        e = zz + M.h[:X]*x #+ M.h[:Z]*z
        mag += S[i, j, m]

        # Z += z
        ZZ += zz
        X += x
        E += e

        for p in 1:P
            p̃ = 2P - p + 2
            if i < L && j < L
                O.D.tmp_∂E[:α][p] += S[i, j, p] * S[i+1, j, p] + S[i, j, p̃] * S[i+1, j, p̃]
                O.D.tmp_∂E[:α][p] += S[i, j, p] * S[i, j+1, p] + S[i, j, p̃] * S[i, j+1, p̃]
            elseif j == L && i < L
                O.D.tmp_∂E[:α][p] += S[i, j, p] * S[i+1, j, p] + S[i, j, p̃] * S[i+1, j, p̃]
            elseif i == L && j < L
                O.D.tmp_∂E[:α][p] += S[i, j, p] * S[i, j+1, p] + S[i, j, p̃] * S[i, j+1, p̃]
            end
            # O.D.tmp_∂E[:γ][p] += S[i, j, p] + S[i, j+1, p̃]
            O.D.tmp_∂E[:β][p] += (S[i, j, p] * S[i, j, p+1] + S[i, j, p̃] * S[i, j, p̃-1])* M.∂Jt[p]
            if p == P
                O.D.tmp_X_∂X[p] += -2S[i, j, m] * S[i, j, m+1] * M.∂Jt[m] * M.h[:X] * x
            end
        end
    end

    return E, ZZ, X, Z, mag
end

function update_thermo!(O::Thermo, c::Array{Int, 3}, M::Model, N::Int, Grad::Gradient)
    P = length(O.D.∂E[:α])
    E, ZZ, X, Z, mag = compute_obs(O, c, M)
    O.E.ZZ += ZZ / N
    # O.E.Z += Z / N
    O.E.X += X / N

    O.E.ZZ² += ZZ^2 / N
    O.E.X² += X^2 / N
    # O.E.Z² += Z^2 / N

    O.E.ZZ_X += ZZ*X / N
    # O.E.ZZ_Z += ZZ*Z / N
    # O.E.X_Z += X*Z / N

    O.E.E += E / N

    O.mag += abs(mag) / N

    for p in 1:P
        for var in [:α, :β]#, :γ]
            O.D.∂E[var][p] += O.D.tmp_∂E[var][p] ./ N
            O.D.E_∂E[var][p] += E * O.D.tmp_∂E[var][p] ./ N
        end
        O.D.E_∂E[:β][p] += O.D.tmp_X_∂X[p] ./ N
    end
end

function update_grad!(grad::Gradient, O::Thermo, M::Model, N::Int)
    N = size(M)[1] * size(M)[2]
    for var in [:α, :β]#, :γ]
        grad.∇E[var] .= (-O.E.E * O.D.∂E[var] .+ O.D.E_∂E[var]) / N
    end
    grad.∇E[:total] .= vec(hcat(grad.∇E[:α], grad.∇E[:β], grad.∇E[:γ])')
end


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

using Markdown
using InteractiveUtils

# ╔═╡ 457e4d2c-0225-11ec-1b21-7b47b5339342
using JuMP, Ipopt

# ╔═╡ 7f3a1eb1-63a7-4062-9dd6-e88e31e961da
model = Model(Ipopt.Optimizer)

# ╔═╡ af9691bd-e434-48e5-bb53-7b63b0d152a9
begin
  @variable(model, x)
  @variable(model, y)
end

# ╔═╡ c199a311-3874-4f2c-a6e9-270cbc40cedd
@NLobjective(model, Min, (1 - x)^2 + 100 * (y - x^2)^2)

# ╔═╡ 8b62d21e-9774-4c5c-9667-c6379e78733c
optimize!(model)

# ╔═╡ 5963c706-4143-443f-aaf0-d61001bfce04
value(x)

# ╔═╡ 2a365ed7-5e15-4b96-851f-2fa556941b51
value(y)

# ╔═╡ 00000000-0000-0000-0000-000000000001
PLUTO_PROJECT_TOML_CONTENTS = """
[deps]
Ipopt = "b6b21f68-93f8-5de0-b562-5493be1d77c9"
JuMP = "4076af6c-e467-56ae-b986-b466b2749572"

[compat]
Ipopt = "~0.7.0"
JuMP = "~0.21.9"
"""

# ╔═╡ 00000000-0000-0000-0000-000000000002
PLUTO_MANIFEST_TOML_CONTENTS = """
# This file is machine-generated - editing it directly is not advised

[[ASL_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "370cafc70604b2522f2c7cf9915ebcd17b4cd38b"
uuid = "ae81ac8f-d209-56e5-92de-9978fef736f9"
version = "0.1.2+0"

[[ArgTools]]
uuid = "0dad84c5-d112-42e6-8d28-ef12dabb789f"

[[Artifacts]]
uuid = "56f22d72-fd6d-98f1-02f0-08ddc0907c33"

[[Base64]]
uuid = "2a0f44e3-6c83-55bd-87e4-b1978d98bd5f"

[[BenchmarkTools]]
deps = ["JSON", "Logging", "Printf", "Statistics", "UUIDs"]
git-tree-sha1 = "aa3aba5ed8f882ed01b71e09ca2ba0f77f44a99e"
uuid = "6e4b80f9-dd63-53aa-95a3-0cdb28fa8baf"
version = "1.1.3"

[[BinaryProvider]]
deps = ["Libdl", "Logging", "SHA"]
git-tree-sha1 = "ecdec412a9abc8db54c0efc5548c64dfce072058"
uuid = "b99e7846-7c00-51b0-8f62-c81ae34c0232"
version = "0.5.10"

[[Bzip2_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "19a35467a82e236ff51bc17a3a44b69ef35185a2"
uuid = "6e34b625-4abd-537c-b88f-471c36dfa7a0"
version = "1.0.8+0"

[[Calculus]]
deps = ["LinearAlgebra"]
git-tree-sha1 = "f641eb0a4f00c343bbc32346e1217b86f3ce9dad"
uuid = "49dc2e85-a5d0-5ad3-a950-438e2897f1b9"
version = "0.5.1"

[[ChainRulesCore]]
deps = ["Compat", "LinearAlgebra", "SparseArrays"]
git-tree-sha1 = "bdc0937269321858ab2a4f288486cb258b9a0af7"
uuid = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4"
version = "1.3.0"

[[CodecBzip2]]
deps = ["Bzip2_jll", "Libdl", "TranscodingStreams"]
git-tree-sha1 = "2e62a725210ce3c3c2e1a3080190e7ca491f18d7"
uuid = "523fee87-0ab8-5b00-afb7-3ecf72e48cfd"
version = "0.7.2"

[[CodecZlib]]
deps = ["TranscodingStreams", "Zlib_jll"]
git-tree-sha1 = "ded953804d019afa9a3f98981d99b33e3db7b6da"
uuid = "944b1d66-785c-5afd-91f1-9de20f533193"
version = "0.7.0"

[[CommonSubexpressions]]
deps = ["MacroTools", "Test"]
git-tree-sha1 = "7b8a93dba8af7e3b42fecabf646260105ac373f7"
uuid = "bbf7d656-a473-5ed7-a52c-81e309532950"
version = "0.3.0"

[[Compat]]
deps = ["Base64", "Dates", "DelimitedFiles", "Distributed", "InteractiveUtils", "LibGit2", "Libdl", "LinearAlgebra", "Markdown", "Mmap", "Pkg", "Printf", "REPL", "Random", "SHA", "Serialization", "SharedArrays", "Sockets", "SparseArrays", "Statistics", "Test", "UUIDs", "Unicode"]
git-tree-sha1 = "727e463cfebd0c7b999bbf3e9e7e16f254b94193"
uuid = "34da2185-b29b-5c13-b0c7-acf172513d20"
version = "3.34.0"

[[CompilerSupportLibraries_jll]]
deps = ["Artifacts", "Libdl"]
uuid = "e66e0078-7015-5450-92f7-15fbd957f2ae"

[[DataStructures]]
deps = ["Compat", "InteractiveUtils", "OrderedCollections"]
git-tree-sha1 = "7d9d316f04214f7efdbb6398d545446e246eff02"
uuid = "864edb3b-99cc-5e75-8d2d-829cb0a9cfe8"
version = "0.18.10"

[[Dates]]
deps = ["Printf"]
uuid = "ade2ca70-3891-5945-98fb-dc099432e06a"

[[DelimitedFiles]]
deps = ["Mmap"]
uuid = "8bb1440f-4735-579b-a4ab-409b98df4dab"

[[DiffResults]]
deps = ["StaticArrays"]
git-tree-sha1 = "c18e98cba888c6c25d1c3b048e4b3380ca956805"
uuid = "163ba53b-c6d8-5494-b064-1a9d43ac40c5"
version = "1.0.3"

[[DiffRules]]
deps = ["NaNMath", "Random", "SpecialFunctions"]
git-tree-sha1 = "3ed8fa7178a10d1cd0f1ca524f249ba6937490c0"
uuid = "b552c78f-8df3-52c6-915a-8e097449b14b"
version = "1.3.0"

[[Distributed]]
deps = ["Random", "Serialization", "Sockets"]
uuid = "8ba89e20-285c-5b6f-9357-94700520ee1b"

[[DocStringExtensions]]
deps = ["LibGit2"]
git-tree-sha1 = "a32185f5428d3986f47c2ab78b1f216d5e6cc96f"
uuid = "ffbed154-4ef7-542d-bbb7-c09d3a79fcae"
version = "0.8.5"

[[Downloads]]
deps = ["ArgTools", "LibCURL", "NetworkOptions"]
uuid = "f43a241f-c20a-4ad4-852c-f6b1247861c6"

[[ForwardDiff]]
deps = ["CommonSubexpressions", "DiffResults", "DiffRules", "LinearAlgebra", "NaNMath", "Printf", "Random", "SpecialFunctions", "StaticArrays"]
git-tree-sha1 = "b5e930ac60b613ef3406da6d4f42c35d8dc51419"
uuid = "f6369f11-7733-5829-9624-2563aa707210"
version = "0.10.19"

[[HTTP]]
deps = ["Base64", "Dates", "IniFile", "Logging", "MbedTLS", "NetworkOptions", "Sockets", "URIs"]
git-tree-sha1 = "44e3b40da000eab4ccb1aecdc4801c040026aeb5"
uuid = "cd3eb016-35fb-5094-929b-558a96fad6f3"
version = "0.9.13"

[[IniFile]]
deps = ["Test"]
git-tree-sha1 = "098e4d2c533924c921f9f9847274f2ad89e018b8"
uuid = "83e8ac13-25f8-5344-8a64-a9f2b223428f"
version = "0.5.0"

[[InteractiveUtils]]
deps = ["Markdown"]
uuid = "b77e0a4c-d291-57a0-90e8-8db25a27a240"

[[Ipopt]]
deps = ["BinaryProvider", "Ipopt_jll", "Libdl", "LinearAlgebra", "MathOptInterface", "MathProgBase"]
git-tree-sha1 = "380786b4929b8d18d76e909c6b2eca355b7c3bd6"
uuid = "b6b21f68-93f8-5de0-b562-5493be1d77c9"
version = "0.7.0"

[[Ipopt_jll]]
deps = ["ASL_jll", "Artifacts", "CompilerSupportLibraries_jll", "JLLWrappers", "Libdl", "MUMPS_seq_jll", "OpenBLAS32_jll", "Pkg"]
git-tree-sha1 = "82124f27743f2802c23fcb05febc517d0b15d86e"
uuid = "9cc047cb-c261-5740-88fc-0cf96f7bdcc7"
version = "3.13.4+2"

[[IrrationalConstants]]
git-tree-sha1 = "f76424439413893a832026ca355fe273e93bce94"
uuid = "92d709cd-6900-40b7-9082-c6be49f344b6"
version = "0.1.0"

[[JLLWrappers]]
deps = ["Preferences"]
git-tree-sha1 = "642a199af8b68253517b80bd3bfd17eb4e84df6e"
uuid = "692b3bcd-3c85-4b1f-b108-f13ce0eb3210"
version = "1.3.0"

[[JSON]]
deps = ["Dates", "Mmap", "Parsers", "Unicode"]
git-tree-sha1 = "8076680b162ada2a031f707ac7b4953e30667a37"
uuid = "682c06a0-de6a-54ab-a142-c8b1cf79cde6"
version = "0.21.2"

[[JSONSchema]]
deps = ["HTTP", "JSON", "ZipFile"]
git-tree-sha1 = "b84ab8139afde82c7c65ba2b792fe12e01dd7307"
uuid = "7d188eb4-7ad8-530c-ae41-71a32a6d4692"
version = "0.3.3"

[[JuMP]]
deps = ["Calculus", "DataStructures", "ForwardDiff", "JSON", "LinearAlgebra", "MathOptInterface", "MutableArithmetics", "NaNMath", "Printf", "Random", "SparseArrays", "SpecialFunctions", "Statistics"]
git-tree-sha1 = "4f0a771949bbe24bf70c89e8032c107ebe03f6ba"
uuid = "4076af6c-e467-56ae-b986-b466b2749572"
version = "0.21.9"

[[LibCURL]]
deps = ["LibCURL_jll", "MozillaCACerts_jll"]
uuid = "b27032c2-a3e7-50c8-80cd-2d36dbcbfd21"

[[LibCURL_jll]]
deps = ["Artifacts", "LibSSH2_jll", "Libdl", "MbedTLS_jll", "Zlib_jll", "nghttp2_jll"]
uuid = "deac9b47-8bc7-5906-a0fe-35ac56dc84c0"

[[LibGit2]]
deps = ["Base64", "NetworkOptions", "Printf", "SHA"]
uuid = "76f85450-5226-5b5a-8eaa-529ad045b433"

[[LibSSH2_jll]]
deps = ["Artifacts", "Libdl", "MbedTLS_jll"]
uuid = "29816b5a-b9ab-546f-933c-edad1886dfa8"

[[Libdl]]
uuid = "8f399da3-3557-5675-b5ff-fb832c97cbdb"

[[LinearAlgebra]]
deps = ["Libdl"]
uuid = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"

[[LogExpFunctions]]
deps = ["DocStringExtensions", "IrrationalConstants", "LinearAlgebra"]
git-tree-sha1 = "3d682c07e6dd250ed082f883dc88aee7996bf2cc"
uuid = "2ab3a3ac-af41-5b50-aa03-7779005ae688"
version = "0.3.0"

[[Logging]]
uuid = "56ddb016-857b-54e1-b83d-db4d58db5568"

[[METIS_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "2dc1a9fc87e57e32b1fc186db78811157b30c118"
uuid = "d00139f3-1899-568f-a2f0-47f597d42d70"
version = "5.1.0+5"

[[MUMPS_seq_jll]]
deps = ["Artifacts", "CompilerSupportLibraries_jll", "JLLWrappers", "Libdl", "METIS_jll", "OpenBLAS32_jll", "Pkg"]
git-tree-sha1 = "1a11a84b2af5feb5a62a820574804056cdc59c39"
uuid = "d7ed1dd3-d0ae-5e8e-bfb4-87a502085b8d"
version = "5.2.1+4"

[[MacroTools]]
deps = ["Markdown", "Random"]
git-tree-sha1 = "0fb723cd8c45858c22169b2e42269e53271a6df7"
uuid = "1914dd2f-81c6-5fcd-8719-6d5c9610ff09"
version = "0.5.7"

[[Markdown]]
deps = ["Base64"]
uuid = "d6f4376e-aef5-505a-96c1-9c027394607a"

[[MathOptInterface]]
deps = ["BenchmarkTools", "CodecBzip2", "CodecZlib", "JSON", "JSONSchema", "LinearAlgebra", "MutableArithmetics", "OrderedCollections", "SparseArrays", "Test", "Unicode"]
git-tree-sha1 = "575644e3c05b258250bb599e57cf73bbf1062901"
uuid = "b8f27783-ece8-5eb3-8dc8-9495eed66fee"
version = "0.9.22"

[[MathProgBase]]
deps = ["LinearAlgebra", "SparseArrays"]
git-tree-sha1 = "9abbe463a1e9fc507f12a69e7f29346c2cdc472c"
uuid = "fdba3010-5040-5b88-9595-932c9decdf73"
version = "0.7.8"

[[MbedTLS]]
deps = ["Dates", "MbedTLS_jll", "Random", "Sockets"]
git-tree-sha1 = "1c38e51c3d08ef2278062ebceade0e46cefc96fe"
uuid = "739be429-bea8-5141-9913-cc70e7f3736d"
version = "1.0.3"

[[MbedTLS_jll]]
deps = ["Artifacts", "Libdl"]
uuid = "c8ffd9c3-330d-5841-b78e-0817d7145fa1"

[[Mmap]]
uuid = "a63ad114-7e13-5084-954f-fe012c677804"

[[MozillaCACerts_jll]]
uuid = "14a3606d-f60d-562e-9121-12d972cd8159"

[[MutableArithmetics]]
deps = ["LinearAlgebra", "SparseArrays", "Test"]
git-tree-sha1 = "3927848ccebcc165952dc0d9ac9aa274a87bfe01"
uuid = "d8a4904e-b15c-11e9-3269-09a3773c0cb0"
version = "0.2.20"

[[NaNMath]]
git-tree-sha1 = "bfe47e760d60b82b66b61d2d44128b62e3a369fb"
uuid = "77ba4419-2d1f-58cd-9bb1-8ffee604a2e3"
version = "0.3.5"

[[NetworkOptions]]
uuid = "ca575930-c2e3-43a9-ace4-1e988b2c1908"

[[OpenBLAS32_jll]]
deps = ["Artifacts", "CompilerSupportLibraries_jll", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "ba4a8f683303c9082e84afba96f25af3c7fb2436"
uuid = "656ef2d0-ae68-5445-9ca0-591084a874a2"
version = "0.3.12+1"

[[OpenSpecFun_jll]]
deps = ["Artifacts", "CompilerSupportLibraries_jll", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "13652491f6856acfd2db29360e1bbcd4565d04f1"
uuid = "efe28fd5-8261-553b-a9e1-b2916fc3738e"
version = "0.5.5+0"

[[OrderedCollections]]
git-tree-sha1 = "85f8e6578bf1f9ee0d11e7bb1b1456435479d47c"
uuid = "bac558e1-5e72-5ebc-8fee-abe8a469f55d"
version = "1.4.1"

[[Parsers]]
deps = ["Dates"]
git-tree-sha1 = "438d35d2d95ae2c5e8780b330592b6de8494e779"
uuid = "69de0a69-1ddd-5017-9359-2bf0b02dc9f0"
version = "2.0.3"

[[Pkg]]
deps = ["Artifacts", "Dates", "Downloads", "LibGit2", "Libdl", "Logging", "Markdown", "Printf", "REPL", "Random", "SHA", "Serialization", "TOML", "Tar", "UUIDs", "p7zip_jll"]
uuid = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f"

[[Preferences]]
deps = ["TOML"]
git-tree-sha1 = "00cfd92944ca9c760982747e9a1d0d5d86ab1e5a"
uuid = "21216c6a-2e73-6563-6e65-726566657250"
version = "1.2.2"

[[Printf]]
deps = ["Unicode"]
uuid = "de0858da-6303-5e67-8744-51eddeeeb8d7"

[[REPL]]
deps = ["InteractiveUtils", "Markdown", "Sockets", "Unicode"]
uuid = "3fa0cd96-eef1-5676-8a61-b3b8758bbffb"

[[Random]]
deps = ["Serialization"]
uuid = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"

[[SHA]]
uuid = "ea8e919c-243c-51af-8825-aaa63cd721ce"

[[Serialization]]
uuid = "9e88b42a-f829-5b0c-bbe9-9e923198166b"

[[SharedArrays]]
deps = ["Distributed", "Mmap", "Random", "Serialization"]
uuid = "1a1011a3-84de-559e-8e89-a11a2f7dc383"

[[Sockets]]
uuid = "6462fe0b-24de-5631-8697-dd941f90decc"

[[SparseArrays]]
deps = ["LinearAlgebra", "Random"]
uuid = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"

[[SpecialFunctions]]
deps = ["ChainRulesCore", "LogExpFunctions", "OpenSpecFun_jll"]
git-tree-sha1 = "a322a9493e49c5f3a10b50df3aedaf1cdb3244b7"
uuid = "276daf66-3868-5448-9aa4-cd146d93841b"
version = "1.6.1"

[[StaticArrays]]
deps = ["LinearAlgebra", "Random", "Statistics"]
git-tree-sha1 = "3240808c6d463ac46f1c1cd7638375cd22abbccb"
uuid = "90137ffa-7385-5640-81b9-e52037218182"
version = "1.2.12"

[[Statistics]]
deps = ["LinearAlgebra", "SparseArrays"]
uuid = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"

[[TOML]]
deps = ["Dates"]
uuid = "fa267f1f-6049-4f14-aa54-33bafae1ed76"

[[Tar]]
deps = ["ArgTools", "SHA"]
uuid = "a4e569a6-e804-4fa4-b0f3-eef7a1d5b13e"

[[Test]]
deps = ["InteractiveUtils", "Logging", "Random", "Serialization"]
uuid = "8dfed614-e22c-5e08-85e1-65c5234f0b40"

[[TranscodingStreams]]
deps = ["Random", "Test"]
git-tree-sha1 = "216b95ea110b5972db65aa90f88d8d89dcb8851c"
uuid = "3bb67fe8-82b1-5028-8e26-92a6c54297fa"
version = "0.9.6"

[[URIs]]
git-tree-sha1 = "97bbe755a53fe859669cd907f2d96aee8d2c1355"
uuid = "5c2747f8-b7ea-4ff2-ba2e-563bfd36b1d4"
version = "1.3.0"

[[UUIDs]]
deps = ["Random", "SHA"]
uuid = "cf7118a7-6976-5b1a-9a39-7adc72f591a4"

[[Unicode]]
uuid = "4ec0a83e-493e-50e2-b9ac-8f72acf5a8f5"

[[ZipFile]]
deps = ["Libdl", "Printf", "Zlib_jll"]
git-tree-sha1 = "c3a5637e27e914a7a445b8d0ad063d701931e9f7"
uuid = "a5390f91-8eb1-5f08-bee0-b1d1ffed6cea"
version = "0.9.3"

[[Zlib_jll]]
deps = ["Libdl"]
uuid = "83775a58-1f1d-513f-b197-d71354ab007a"

[[nghttp2_jll]]
deps = ["Artifacts", "Libdl"]
uuid = "8e850ede-7688-5339-a07c-302acd2aaf8d"

[[p7zip_jll]]
deps = ["Artifacts", "Libdl"]
uuid = "3f19e933-33d8-53b3-aaab-bd5110c3b7a0"
"""

# ╔═╡ Cell order:
# ╠═457e4d2c-0225-11ec-1b21-7b47b5339342
# ╠═7f3a1eb1-63a7-4062-9dd6-e88e31e961da
# ╠═af9691bd-e434-48e5-bb53-7b63b0d152a9
# ╠═c199a311-3874-4f2c-a6e9-270cbc40cedd
# ╠═8b62d21e-9774-4c5c-9667-c6379e78733c
# ╠═5963c706-4143-443f-aaf0-d61001bfce04
# ╠═2a365ed7-5e15-4b96-851f-2fa556941b51
# ╟─00000000-0000-0000-0000-000000000001
# ╟─00000000-0000-0000-0000-000000000002


text:
using
    GPUAcceleratedTracking,
    Test,
    Tracking,
    GNSSSignals,
    CUDA,
    StructArrays

import Tracking: Hz, ms

if CUDA.functional()
    include("algorithms.jl")
    include("reduction.jl")
end

max_stars_repo_name: mfkiwl/GPUAcceleratedTracking

max_issues_repo_name: mfkiwl/GPUAcceleratedTracking

max_forks_repo_name: mfkiwl/GPUAcceleratedTracking


text:
function Vth_EKV(Vfb,Na,Cox,T,phif)
    GAMMA = γC(Na,Cox)
    PHI = 2*phif
    VTO = Vfb+PHI+GAMMA*sqrt(PHI)
end

function Vth_EKV(mos::MOSStructure)
    Vth_EKV(VFB(mos),abs(C(mos)),Cox(mos),Temperature(mos),ϕb(mos))
end

function Vp_EKV(Vg,Vfb,Na,Cox,T,phif,sigma =0.0 )
    GAMMA = γC(Na,Cox)
    PHI = 2*phif
    VTO = Vfb+PHI+GAMMA*sqrt(PHI)
    VGprime = Vg - VTO + PHI + GAMMA * sqrt(PHI)
    VP = VGprime - PHI - GAMMA* (sqrt(VGprime+(GAMMA/2.0)*(GAMMA/2.0))-(GAMMA/2.0));
    return VP
end

function Vp_EKV(Vg,mos::MOSStructure )
    Vp_EKV(Vg,VFB(mos),abs(C(mos)),Cox(mos),Temperature(mos),ϕb(mos))
end

function n_EKV(Vg,Vfb,Na,Cox,T,phif,sigma =0.0 )
    GAMMA = γC(Na,Cox)
    PHI = 2*phif
    vt = kb*T
    VP = Vp_EKV(Vg,Vfb,Na,Cox,T,phif)
    n = 1.0 + GAMMA / (2.0*sqrt(PHI + VP + 4.0*vt));
    return n
end

function n_EKV(Vg,mos::MOSStructure )
    n_EKV(Vg,VFB(mos),abs(C(mos)),Cox(mos),Temperature(mos),ϕb(mos))
end

function G_EKV(Na,Cox,T)
    γC(Na,Cox)/sqrt(kb*T)
end

function G_EKV(mos::MOSStructure)
    γC(abs(C(mos)),Cox(mos))/sqrt(kb*Temperature(mos))
end

function q_EKV(Vg,V,Vfb,Na,Cox,T,phif,sigma =0.0 )
    VP = Vp_EKV(Vg,Vfb,Na,Cox,T,phif)
    x=(VP-V)/kb/T
    iff = (log(1.0+exp( x /2.0)))*(log(1.0+exp( x /2.0)));
    return iff
end

function q_EKV(Vg,V,mos::MOSStructure)
   q_EKV(Vg,V,VFB(mos),abs(C(mos)),Cox(mos),Temperature(mos),ϕb(mos))
end

function id_EKV(Vg,Vd,Vs,Vfb,Na,Cox,T,phif,sigma =0.0 )
    pt = kb*T
    Ispec = 2*n_EKV(Vg,Vfb,Na,Cox,T,phif)*pt*pt*Cox
    return Ispec*(q_EKV(Vg,Vs,Vfb,Na,Cox,T,phif,sigma)-q_EKV(Vg,Vd,Vfb,Na,Cox,T,phif,sigma))

end

function Id_EKV(Vg,Vd,Vs,Vfb,Na,Cox,T,phif,sigma =0.0,Kp=5e-2 )
    pt = kb*T
    Ispec = 2*n_EKV(Vg,Vfb,Na,Cox,T,phif,sigma)*pt*pt*Cox
    Kp = Kp*(T/300)^(-1.5)
    return Kp*Ispec*(q_EKV(Vg,Vs,Vfb,Na,Cox,T,phif,sigma)-q_EKV(Vg,Vd,Vfb,Na,Cox,T,phif,sigma))

end

mutable struct EKV26Model <: TransistorModel
    mos::MOSStructure
    W::Number
    L::Number
    μ_0
end

id(vg,vd,vs,m::EKV26Model) = id_EKV(vg,vd,vs,VFB(m.mos),abs(C(m.mos)),Cox(m.mos),Temperature(m.mos),ϕb(m.mos))
Id(vg,vd,vs,m::EKV26Model,sigma=0.0,Kp=5e-2) = Id_EKV(vg,vd,vs,VFB(m.mos),abs(C(m.mos)),Cox(m.mos),Temperature(m.mos),ϕb(m.mos),sigma,Kp)
Vp(Vg,m::EKV26Model) = Vp_EKV(Vg,m.mos)


text:
module Badges

using Genie, Stipple, StippleUI, StippleUI.API
import Genie.Renderer.Html: HTMLString, normal_element

export badge, Badge

function __init__()
  Genie.Renderer.Html.register_normal_element("q__badge", context = Genie.Renderer.Html)
end

"""
floating::Bool = false
transparent::Bool = false
multiline::Bool = true
label::String = "Label"
align::String = "top" # middle, bottom
outline::Bool = true
color::String = "blue"
textcolor::String = "teal-10"
"""
function badge( fieldname::Union{Symbol,Nothing} = nothing,
                args...;
                wrap::Function = StippleUI.DEFAULT_WRAPPER,
                kwargs...) where {T<:Stipple.ReactiveModel}
  wrap() do
    Genie.Renderer.Html.q__badge(args...;
            attributes(
              [:fieldname => fieldname, kwargs...],
              StippleUI.API.ATTRIBUTES_MAPPINGS
            )...)
  end
end

mutable struct Badge
  fieldname
  args
  wrap
  kwargs

  Badge(fieldname::Union{Symbol,Nothing} = nothing,
        args...;
        wrap::Function = StippleUI.DEFAULT_WRAPPER,
        kwargs...) = new(fieldname, args, wrap, kwargs)
end

Base.string(b::Badge) = badge(b.fieldname, b.args...; b.wrap, b.kwargs...)

end


text:
#= single parameter error-free transformations =#

# 'y' must be negated to get the right result
@inline function eftRecip{T<:AbstractFloat}(a::T)
     x = one(T)/a
     y = -(fma(x,a,-1.0)/a)
     x,y
end

#=
   While not strictly an error-free transformation it is quite reliable and recommended for use.
   Augmented precision square roots, 2-D norms and discussion on correctly reounding sqrt(x^2 + y^2)
   by Nicolas Brisebarre, Mioara Joldes, Erik Martin-Dorel, Hean-Michel Muller, Peter Kornerup
=#

@inline function eftSqrt(a::AbstractFloat)
     x = sqrt(a)
     t = fma(x,-x,a)
     y = t / (x*2.0)
     x,y
end

@inline function eftRecipSqrt(a::Float64)
     r = 1.0/a
     x = sqrt(r)
     t = fma(x,-x,r)
     y = t / (x*2.0)
     x,y
end

@inline function eftSquare(a::Float64)
    x = a * a
    y = fma(a, a, -x)
    x,y
end

@inline function eftCube(a::Float64)
    p = a*a; e = fma(a, a, -p)
    x = p*a; p = fma(p, a, -x)
    y = e*a
    x,y
end

#= two parameter error-free transformations =#

@inline function eftSum2{T<:AbstractFloat}(a::T, b::T)
  x = a + b
  t = x - a
  y = (a - (x - t)) + (b - t)
  x,y
end

@inline function eftSum2inOrder{T<:AbstractFloat}(a::T, b::T)
  x = a + b
  y = b - (x - a)
  x,y
end

@inline function eftDiff2{T<:AbstractFloat}(a::T, b::T)
  x = a - b
  t = x - a
  y = (a - (x - t)) - (b + t)
  x,y
end

@inline function eftDiff2inOrder{T<:AbstractFloat}(a::T, b::T)
  x = a - b
  y = (a - x) - b
  x,y
end

@inline function eftProd2{T<:AbstractFloat}(a::T, b::T)
    x = a * b
    y = fma(a, b, -x)
    x,y
end

# 'y' must be negated to get the right result
@inline function eftDiv2{T<:AbstractFloat}(a::T,b::T)
     x = a/b
     y = -(fma(x,b,-a)/b)
     x,y
end

#= three parameter error-free transformations =#

function eftSum3{T<:Float64}(a::T,b::T,c::T)
    s,t = eftSum2(b, c)
    x,u = eftSum2(a, s)
    y,z = eftSum2(u, t)
    x,y = eftSum2inOrder(x, y)
    x,y,z
end

function eftSum3inOrder{T<:Float64}(a::T,b::T,c::T)
    s,t = eftSum2inOrder(b, c)
    x,u = eftSum2inOrder(a, s)
    y,z = eftSum2inOrder(u, t)
    x,y = eftSum2inOrder(x, y)
    x,y,z
end

function eftProd3{T<:Float64}(a::T, b::T, c::T)
    p,e = eftProd2(a,b)
    x,p = eftProd2(p,c)
    y,z = eftProd2(e,c)
    x,y,z
end

function eftFMA{T<:Float64}(a::T, b::T, c::T)
    x = fma(a,b,c)
    u1,u2 = eftProd2(a,b)
    a1,a2 = eftSum2(u2,c)
    b1,b2 = eftSum2(u1,a1)
    g = (b1-x)+b2
    y,z = eftSum2inOrder(g,a2)
    x,y,z
end

function eftFMS{T<:Float64}(a::T, b::T, c::T)
    x = fma(a,b,c)
    u1,u2 = eftProd2(a,b)
    a1,a2 = eftDiff2(u2,c)
    b1,b2 = eftSum2(u1,a1)
    g = (b1-x)+b2
    y,z = eftSum2inOrder(g,a2)
    x,y,z
end


text:
module Test_bouncingBall_as_dae

#=
DAE solver using IDA and events: solution sol seem to store wrong values

https://github.com/JuliaDiffEq/DifferentialEquations.jl/issues/552
=#

using DifferentialEquations, PyPlot, Sundials

function derivative(du,u,p,t)
    du[1] = u[2]
    du[2] = -9.81
end

function f(res,du,u,p,t)
  derivative(res,u,p,t)
  res[1] = du[1] - res[1]
  res[2] = du[2] - res[2]
end

function condition(u,t,integrator) # Event when event_f(u,t) == 0
 z = u[1]  + 1e-12
 println("condition called at time = ", t, ", z = ", z)
 return z
end

function affect_neg(integrator)
  println("time = ", integrator.t, ", h = ", integrator.u[1])
  integrator.u[2] = -0.8*integrator.u[2]
  derivative(integrator.du, integrator.u, integrator.p, integrator.t)
end

cb = ContinuousCallback(condition,nothing,affect_neg! = affect_neg)

u0 = [1.0,0.0]
du0= [0.0,-9.81]
tspan = (0.0,1.5)
prob = DAEProblem(f,du0,u0,tspan)
sol = solve(prob,IDA(),saveat=0.01,callback=cb)

figure(1)
clf()
plot(sol.t, getindex.(sol.u, 1), label="h")

end


text:
# From Baltensperger and Trummer (eqns 3 - 5) we have the following.
function D1(r::T) where {T <: BRInterp}
    n = length(r.x)
    D = zeros(eltype(r.x), n, n)
    @inbounds for j ∈ 1:n
        for i ∈ 1:n
            if i != j
                D[i, j] = r.w[j] / (r.w[i] * (r.x[i] - r.x[j]))
            end
        end
        D[j, j] = - r.x[j] / (2*(one(eltype(r.x)) - r.x[j]^2))
    end
    D[1, 1] = (2*(n-1)^2 + 1) / 6
    D[end, end] = - D[1,1]
    return D
end

function differ(r::T) where {T <: FHInterp}
    n = length(r.x)

    # First the polynomial contribution
    D = D1(r)
    df = D * r.f

    # Now the rational "correction"
    dd2 = divided_difference(r.x ,r.f) ^ 2
    t = zeros(eltype(r.f), n)
    @inbounds for j ∈ 1:n
        t .= r.x[j] .- r.x
        idx = [1:j-1; j+1;n;]
        @views xv = r.x[idx]
        @views fv = r.f[idx]
        t = t ./ divided_difference(xv , fv)
        t[j] = one(eltype(r.x))
        df[j] += dd2 * prod(t)
    end

    dr = FHInterp(r.x, df, r.w, r.order)
end

function differ(r::T) where {T <: AAAapprox}
    D = D1(r)
    df = D * r.f
    dr = AAAapprox(r.x, df, r.w, zeros(eltype(r.x), 1))
end

function divided_difference(x, y)
    n = length(x)
    s = zero(eltype(x))
    t = zeros(eltype(x), n)
    @inbounds for j ∈ 1:n
        t .= x[j] .- x
        t[j] = one(eltype(x))
        s += y[j] / prod(t)
    end
    return s
end


text:
#!/usr/bin/env julia

using Luxor, Colors

function test_axes(x, y, rot, w)
    gsave()
    translate(x, y)
    rotate(rot)
    box(0, 0, w, w, :clip)
    axes()
    clipreset()
    grestore()
end

fname = "/tmp/axes-test.pdf"
width, height = 2000, 2000
Drawing(width, height, fname)
origin()
background("ivory")
pagetiles = Tiler(width, height, 5, 5, margin=50)
for (pos, n) in pagetiles
  test_axes(pos.x, pos.y, rand() * 2pi, pagetiles.tilewidth)
end

finish()
println("finished test: output in $(fname)")


text:
using Flux
include("hparams.jl")

function Discriminator()
	return Chain(
			x -> reshape(x,4,:),
			Dense(4,128),
			x -> leakyrelu.(x,0.3f0),
			Dropout(0.2),
			Dense(128, 64),
			x -> leakyrelu.(x, 0.2f0),
			Dropout(0.2),
			x -> reshape(x,2,2,16,:),
			Conv((4, 4), 16 => 8; stride = 1, pad = 1),
			x->leakyrelu.(x,0.2f0),
			Dropout(0.2),
			Conv((3,3), 8=>4; stride = 1, pad = 1),
			x->leakyrelu.(x,0.2f0),
			x->reshape(x,3*4,:),
			Dense(3*4, 1),
			x->leakyrelu.(x,0.2f0),
			)
end

function Generator(hparams)
     return Chain(
			x -> reshape(x,4,4,12,:),
			Conv((2,2), 12=>hparams.latent_dim;stride=1,pad=1),
			x -> leakyrelu.(x,0.2f0),
			Dropout(0.2),
			x -> reshape(x,hparams.latent_dim,:),
            Dense(hparams.latent_dim, 3*3*64),
			x -> leakyrelu.(x,0.3f0),
            Dropout(0.3),
            BatchNorm(5 * 5 * 64, relu),
            x->reshape(x, 5, 5, 64, :),
            ConvTranspose((3, 3), 64 => 16; stride = 1, pad = 2),
            BatchNorm(16, relu),
            ConvTranspose((3, 3), 16 => 8; stride = 2, pad = 1),
            BatchNorm(8, relu),
            ConvTranspose((3, 3), 8 => 3, tanh; stride = 2, pad = 1),
            x->transpose.(x)
			)
end

function Observer()
     return Chain(
             x -> reshape(x,4,:),
             Dense(4,256),
             x -> leakyrelu.(x,0.15f0),
             Dropout(0.25),
             BatchNorm(256, relu),
             x -> reshape(x, 4,4,16,:),
             Conv((3, 3), 16 => 64; stride = 1, pad = 1),
             x -> leakyrelu.(x,0.2f0),
             Dropout(0.2),
			 x -> reshape(x,64,:),
             Dense(64, 1),
             x -> leakyrelu.(x,0.2f0))
end


text:
module CGP

# This defines what the import module of 'CGP is'

using Logging
using PaddedViews
using Distributions

Logging.configure(level=INFO)

# Various includes we wish to have in this module

include("config.jl")
include("chromosome.jl")
include("distance.jl")
include("mutation.jl")
include("crossover.jl")
include("chromosomes/cgp.jl")
include("chromosomes/pcgp.jl")
include("logging.jl")
include("EAs/oneplus.jl")
include("EAs/cgpneat.jl")
include("EAs/ga.jl")
include("EAs/cmaes.jl")

# The series of evolutionary strategies to be included and used
EAs = [oneplus, GA]

#  The types of chromosomes, normal CGP and positional CGP
chromosomes = [CGPChromo, PCGPChromo]

# The types of mutation strategies that are included in the package
mutations = [:gene_mutate, :mixed_node_mutate, :mixed_subtree_mutate, :adaptive_node_mutate,
             :adaptive_subtree_mutate]

# The types of crossover strategies included in the package
crossovers = [:single_point_crossover, :random_node_crossover, :aligned_node_crossover,
              :proportional_crossover, :output_graph_crossover, :subgraph_crossover]

# The distance measures included in the package
distances = [:genetic_distance, :positional_distance, :constant_functional_distance,
             :random_functional_distance, :active_distance]

end


text:
# This file contains code that was formerly a part of Julia.
# License is MIT: https://julialang.org/license

module KahanSummation

export sum_kbn, cumsum_kbn

"""
    cumsum_kbn(A; dims::Integer)

Cumulative sum along a dimension, using the Kahan-Babuska-Neumaier compensated summation
algorithm for additional accuracy.
"""
function cumsum_kbn(A::AbstractArray{T}; dims::Integer) where T<:AbstractFloat
    axis_size = size(A, dims)
    axis_stride = 1
    for i = 1:dims-1
        axis_stride *= size(A, i)
    end
    axis_size <= 1 && return A
    B = similar(A)
    C = similar(A)
    for i = 1:length(A)
        if div(i-1, axis_stride) % axis_size == 0
            B[i] = A[i]
            C[i] = zero(T)
        else
            s = B[i-axis_stride]
            Ai = A[i]
            B[i] = t = s + Ai
            if abs(s) >= abs(Ai)
                C[i] = C[i-axis_stride] + ((s-t) + Ai)
            else
                C[i] = C[i-axis_stride] + ((Ai-t) + s)
            end
        end
    end
    return B + C
end

function cumsum_kbn(v::AbstractVector{T}) where T<:AbstractFloat
    r = similar(v)
    isempty(v) && return r
    inds = axes(v, 1)
    i1 = first(inds)
    s = r[i1] = v[i1]
    c = zero(T)
    for i = i1+1:last(inds)
        vi = v[i]
        t = s + vi
        if abs(s) >= abs(vi)
            c += ((s-t) + vi)
        else
            c += ((vi-t) + s)
        end
        s = t
        r[i] = s+c
    end
    return r
end

"""
    sum_kbn(A)

Return the sum of all elements of `A`, using the Kahan-Babuska-Neumaier compensated
summation algorithm for additional accuracy.
"""
function sum_kbn(A)
    T = Base.@default_eltype(A)
    c = Base.reduce_empty(+, T)
    it = iterate(A)
    it === nothing && return c
    Ai, i = it
    s = Ai - c
    while (it = iterate(A, i)) !== nothing
        Ai, i = it
        t = s + Ai
        if abs(s) >= abs(Ai)
            c -= ((s-t) + Ai)
        else
            c -= ((Ai-t) + s)
        end
        s = t
    end
    s - c
end

### Deprecations

Base.@deprecate cumsum_kbn(A, axis) cumsum_kbn(A, dims=axis)

end # module


text:
using Test

@testset "FEMfunctions.jl" begin
  include("test_onedElliptic.jl")

  include("test_twodElliptic.jl")

  include("test_TriMesh.jl")

  include("test_twodAdvectionDiffusion.jl")
end


text:
# Effective Rate Factor
ERF(::ContinuousCompounding, r::Float64, t::Float64) = t == 0.0 ? 1.0 : exp(r*t)
ERF(::SimpleCompounding, r::Float64, t::Float64) = t == 0.0 ? 1.0 : 1.0 + r*t
ERF(::ExponentialCompounding, r::Float64, t::Float64) = t == 0.0 ? 1.0 : (1.0+r)^t
ERF(ct::CompoundingType, dcc::DayCountConvention, r::Float64, date_start::Date, date_end::Date) = ERF(ct, r, yearfraction(dcc, date_start, date_end))
ERF(curve::AbstractIRCurve, maturity::Date) = ERF(curve_get_compounding(curve), curve_get_daycount(curve), zero_rate(curve, maturity), curve_get_date(curve), maturity)
ERF(curve::AbstractIRCurve, forward_date::Date, maturity::Date) = ERF(curve, maturity) / ERF(curve, forward_date)

ERF_to_rate(::ContinuousCompounding, ERF::Float64, t::Float64) = log(ERF) / t
ERF_to_rate(::SimpleCompounding, ERF::Float64, t::Float64) = (ERF-1.0) / t
ERF_to_rate(::ExponentialCompounding, ERF::Float64, t::Float64) = ERF^(1.0/t) - 1.0
ERF_to_rate(curve::AbstractIRCurve, ERF::Float64, t::Float64) = ERF_to_rate(curve_get_compounding(curve), ERF, t)

discountfactor_to_rate(c::CompoundingType, _discountfactor_::Float64, t::Float64) = ERF_to_rate(c, 1.0 / _discountfactor_, t)

function discountfactor_to_rate(c::CompoundingType, _discountfactor_vec_::Vector{Float64}, t_vec::Vector{Float64})
    l = length(_discountfactor_vec_)
    @assert l == length(t_vec) "_discountfactor_vec_ and t_vec must have the same length. ($l != $(length(t_vec)))"

    result = Vector{Float64}(undef, l)
    for i in 1:l
        @inbounds result[i] = discountfactor_to_rate(c, _discountfactor_vec_[i], t_vec[i])
    end

    return result
end

# Effective Rate = [Effective Rate Factor] - 1
ER(c::CompoundingType, r::Float64, t::Float64) = ERF(c, r, t) - 1.0
ER(ct::CompoundingType, dcc::DayCountConvention, r::Float64, date_start::Date, date_end::Date) = ER(ct, r, yearfraction(dcc, date_start, date_end))
ER(curve::AbstractIRCurve, maturity::Date) = ERF(curve, maturity) - 1.0
ER(curve::AbstractIRCurve, forward_date::Date, maturity::Date) =  ERF(curve, forward_date, maturity) - 1.0

# [Discount Factor] = 1 / [Effective Rate Factor]
discountfactor(ct::CompoundingType, r::Float64, t::Float64) = 1.0 / ERF(ct, r, t)
discountfactor(ct::CompoundingType, dcc::DayCountConvention, r::Float64, date_start::Date, date_end::Date) = discountfactor(ct, r, yearfraction(dcc, date_start, date_end))
discountfactor(curve::AbstractIRCurve, maturity::Date) = 1.0 / ERF(curve, maturity)

# Optimized vector functions for `ERF` and `discountfactor` functions
for fun in (:ERF, :discountfactor)
    @eval begin
        function ($fun)(curve::AbstractIRCurve, maturity_vec::Vector{Date})
            len = length(maturity_vec)
            _zero_rate_vec_ = zero_rate(curve, maturity_vec)
            result = Vector{Float64}(undef, len)
            for i = 1:len
                @inbounds result[i] = $(fun)(curve_get_compounding(curve), curve_get_daycount(curve), _zero_rate_vec_[i], curve_get_date(curve), maturity_vec[i])
            end
            return result
        end
    end
end


text:
@testset "preprocess" begin

@testset "all_hydrogens" begin
    ethanol = smilestomol("[H]C([H])([H])C([H])([H])O")
    @test issetequal(allhydrogens(ethanol), [1, 3, 4, 6, 7])
    removed = makehydrogensimplicit(ethanol)
    @test nodecount(removed) == 3
    @test edgecount(removed) == 2
end

@testset "trivial_hydrogens" begin
    ethanol = smilestomol("[H]C([H])([H])C([H])([H])O")
    @test issetequal(trivialhydrogens(ethanol), [1, 3, 4, 6, 7])
    removed = makehydrogensimplicit(ethanol, all=false)
    @test nodecount(removed) == 3
    @test edgecount(removed) == 2
end

@testset "hydrogens_explicit" begin
    neop = smilestomol("CC(C)(C)CO")
    added = makehydrogensexplicit(neop)
    @test nodecount(added) == 18
    @test edgecount(added) == 17
end

@testset "largest_component" begin
    thiols = smilestomol("CCCCS.SCCCCC")
    @test issetequal(largestcomponentnodes(thiols), 6:11)
    lc = largestcomponentgraph(thiols)
    @test nodecount(lc) == 6
    @test edgecount(lc) == 5
end

@testset "neutralize_acids" begin
    AcOH = smilestomol("CC(=O)[O-]")
    neutralizeacids!(AcOH)
    @test nodeattr(AcOH, 4).charge == 0
end

@testset "neutralize_oniums" begin
    pyrrolidinium = smilestomol("C1[N+]CCC1")
    neutralizeoniums!(pyrrolidinium)
    @test nodeattr(pyrrolidinium, 2).charge == 0
end

@testset "depolarize" begin
    acetone = smilestomol("C[C+]([O-])C")
    depolarize!(acetone)
    @test nodeattr(acetone, 2).charge == 0
    @test nodeattr(acetone, 3).charge == 0
    @test edgeattr(acetone, 2).order == 2
end

@testset "triplebond_anion" begin
    me_azide = smilestomol("C[N-][N+]#N")
    triplebondanion!(me_azide)
    @test nodeattr(me_azide, 2).charge == 0
    @test nodeattr(me_azide, 4).charge == -1
    @test edgeattr(me_azide, 2).order == 2
    @test edgeattr(me_azide, 3).order == 2
end

end # preprocess


text:
function helper_fn(z, attempt)
	single_data_name = z
	proc_id = myid();

	helper_start_time = time()

	choice_probs = readcsv( string("../tcm_data/synthetic_data_models/", single_data_name, "_model_predict.csv") )
	grand_choice_probs = readcsv( string("../tcm_data/synthetic_data_models/", single_data_name, "_grand_model_predict.csv") )

	train_assortments = convert(Array{Int64}, readcsv( string("../tcm_data/synthetic_data_assortments/", single_data_name, "_assortments.csv") ) )

	N = 10
	grand_assortments = tcm_createGrandAssortments(N)
	time_limit = 1200

	num_train_assortments = size(train_assortments,1)
	num_test_assortments = size(grand_assortments, 1)

	revenues = collect(90:-10:0)

	grand_best_r = maximum( grand_choice_probs * revenues )

	## RANKING-BASED MODEL ##
	numInitializations_LS = 5
	ranking_lambda, orderings, ranking_L1_train_error, ranking_elapsed_time = tcm_estimateRanking_L1(N, train_assortments, choice_probs, numInitializations_LS, time_limit )
	ranking_K = length(ranking_lambda)

	ranking_grand_predict, ranking_grand_predict_by_class = tcm_predictRanking(N, ranking_K, ranking_lambda, orderings, grand_assortments)
	ranking_L1_test_error = sum( abs.( ranking_grand_predict - grand_choice_probs  )  ) / 2^(N-1)

	ranking_best_a = indmax( ranking_grand_predict * revenues)
	ranking_best_r = dot(grand_choice_probs[ranking_best_a,:], revenues)

	ranking_approx_rate = ranking_best_r / grand_best_r

	row_string = string(single_data_name,  ",",
							num_train_assortments, ",",
							num_test_assortments, ",",
							proc_id, ",",
							attempt, ",",
							"ranking", ",",
							ranking_L1_train_error, ",",
							ranking_K, ",",
							ranking_L1_test_error, ",",
							ranking_approx_rate, ",",
							ranking_elapsed_time)
	row_string = string(row_string, "\n");

	@spawnat 1 begin
		print(ranking_outcsvhandle, row_string);
		flush(ranking_outcsvhandle);
	end

	## LC-MNL MODEL ##
	perTol = 1e-6
	K_list = [5]
	numInitializations = 5

	for K in K_list
		LCMNL_u_val, LCMNL_class_probs, LCMNL_loglik, LCMNL_grand_elapsed_time = tcm_estimateMMNL(N, K, train_assortments, choice_probs, numInitializations, time_limit, perTol)

		LCMNL_train_predict, LCMNL_train_predict_by_class = tcm_predictMMNL(N, K, LCMNL_class_probs, LCMNL_u_val, train_assortments)
		LCMNL_L1_train_error = sum( abs.( LCMNL_train_predict - choice_probs  )  )

		LCMNL_grand_predict, LCMNL_grand_predict_by_class = tcm_predictMMNL(N, K, LCMNL_class_probs, LCMNL_u_val, grand_assortments)
		LCMNL_L1_test_error = sum( abs.( LCMNL_grand_predict - grand_choice_probs  )  ) / 2^(N-1)

		LCMNL_best_a = indmax( LCMNL_grand_predict * revenues)
		LCMNL_best_r = dot(grand_choice_probs[LCMNL_best_a,:], revenues)

		LCMNL_approx_rate = LCMNL_best_r / grand_best_r

		row_string = string(single_data_name,  ",",
							num_train_assortments, ",",
							num_test_assortments, ",",
							proc_id, ",",
							attempt, ",",
							"LCMNL", ",",
							LCMNL_loglik, ",",
							K, ",",
							perTol, ",",
							numInitializations, ",",
							LCMNL_L1_train_error, ",",
							LCMNL_L1_test_error, ",",
							LCMNL_approx_rate, ",",
							LCMNL_grand_elapsed_time)
		row_string = string(row_string, "\n");

		@spawnat 1 begin
			print(MNL_outcsvhandle, row_string);
			flush(MNL_outcsvhandle);
		end
	end

	## MNL MODEL ##

	MNL_u_val, MNL_loglik, MNL_elapsed_time = tcm_estimateMNL(N, train_assortments, choice_probs)
	MNL_train_predict = tcm_predictMNL(N, MNL_u_val, train_assortments)
	MNL_L1_train_error = sum( abs.( MNL_train_predict - choice_probs  )  )

	MNL_grand_predict = tcm_predictMNL(N, MNL_u_val, grand_assortments)
	MNL_L1_test_error = sum( abs.( MNL_grand_predict - grand_choice_probs  )  ) / 2^(N-1)

	MNL_best_a = indmax( MNL_grand_predict * revenues)
	MNL_best_r = dot(grand_choice_probs[MNL_best_a,:], revenues)
	MNL_approx_rate = MNL_best_r / grand_best_r

	row_string = string(single_data_name,  ",",
							num_train_assortments, ",",
							num_test_assortments, ",",
							proc_id, ",",
							attempt, ",",
							"MNL", ",",
							MNL_loglik, ",",
							1, ",",
							0, ",",
							1, ",",
							MNL_L1_train_error, ",",
							MNL_L1_test_error, ",",
							MNL_approx_rate, ",",
							MNL_elapsed_time)
	row_string = string(row_string, "\n");

	@spawnat 1 begin
		print(MNL_outcsvhandle, row_string);
		flush(MNL_outcsvhandle);
	end

	## HALO-MNL MODEL ##

	HALOMNL_u_val, HALOMNL_alpha_val, HALOMNL_loglik, HALOMNL_elapsed_time = tcm_estimateHALOMNL(N, train_assortments, choice_probs)
	HALOMNL_train_predict = tcm_predictHALOMNL(N, HALOMNL_u_val, HALOMNL_alpha_val, train_assortments)
	HALOMNL_L1_train_error = sum( abs.( HALOMNL_train_predict - choice_probs  )  )

	HALOMNL_grand_predict = tcm_predictHALOMNL(N, HALOMNL_u_val, HALOMNL_alpha_val, grand_assortments)
	HALOMNL_L1_test_error = sum( abs.( HALOMNL_grand_predict - grand_choice_probs  )  ) / 2^(N-1)

	HALOMNL_best_a = indmax( HALOMNL_grand_predict * revenues)
	HALOMNL_best_r = dot(grand_choice_probs[MNL_best_a,:], revenues)
	HALOMNL_approx_rate = HALOMNL_best_r / grand_best_r

	row_string = string(single_data_name,  ",",
							num_train_assortments, ",",
							num_test_assortments, ",",
							proc_id, ",",
							attempt, ",",
							"HALOMNL", ",",
							HALOMNL_loglik, ",",
							1, ",",
							0, ",",
							1, ",",
							HALOMNL_L1_train_error, ",",
							HALOMNL_L1_test_error, ",",
							HALOMNL_approx_rate, ",",
							HALOMNL_elapsed_time)
	row_string = string(row_string, "\n");

	@spawnat 1 begin
		print(MNL_outcsvhandle, row_string);
		flush(MNL_outcsvhandle);
	end

	## FOREST, COLD-START ##

	depth_list = [4]

	initial_left, initial_right, initial_product, initial_isLeaf = tcm_createBasicForest(N)

	for depth_limit in depth_list

		lambda, forest_left, forest_right, forest_product, forest_isLeaf, forest_L1_train_error, forest_elapsed_time = tcm_estimateForest_L1(N, train_assortments, choice_probs, time_limit, depth_limit,
																												initial_left, initial_right, initial_product, initial_isLeaf)
		K = length(lambda)

		forest_grand_predict, forest_grand_predict_by_class = tcm_predictForest(N, K, lambda, forest_left, forest_right, forest_product, forest_isLeaf, grand_assortments)
		forest_L1_test_error = sum( abs.( forest_grand_predict - grand_choice_probs  )  ) / 2^(N-1)

		forest_best_a = indmax( forest_grand_predict * revenues)
		forest_best_r = dot(grand_choice_probs[forest_best_a,:], revenues)
		forest_approx_rate = forest_best_r / grand_best_r

		row_string = string(single_data_name,  ",",
							num_train_assortments, ",",
							num_test_assortments, ",",
							proc_id, ",",
							attempt, ",",
							"forest-cold", ",",
							depth_limit, ",",
							forest_L1_train_error, ",",
							K, ",",
							forest_L1_test_error, ",",
							forest_approx_rate, ",",
							forest_elapsed_time)
		row_string = string(row_string, "\n");

		@spawnat 1 begin
			print(forest_outcsvhandle, row_string);
			flush(forest_outcsvhandle);
		end
	end

	## FOREST, COLD-START (EMPTY INITIAL FOREST) ##

	depth_list = [4]

	initial_left, initial_right, initial_product, initial_isLeaf = Array{Int64,1}[], Array{Int64,1}[], Array{Int64,1}[], Array{Bool,1}[]; #tcm_createBasicForest(N)

	for depth_limit in depth_list

		lambda, forest_left, forest_right, forest_product, forest_isLeaf, forest_L1_train_error, forest_elapsed_time = tcm_estimateForest_L1(N, train_assortments, choice_probs, time_limit, depth_limit,
																												initial_left, initial_right, initial_product, initial_isLeaf)
		K = length(lambda)

		forest_grand_predict, forest_grand_predict_by_class = tcm_predictForest(N, K, lambda, forest_left, forest_right, forest_product, forest_isLeaf, grand_assortments)
		forest_L1_test_error = sum( abs.( forest_grand_predict - grand_choice_probs  )  ) / 2^(N-1)

		forest_best_a = indmax( forest_grand_predict * revenues)
		forest_best_r = dot(grand_choice_probs[forest_best_a,:], revenues)
		forest_approx_rate = forest_best_r / grand_best_r

		row_string = string(single_data_name,  ",",
							num_train_assortments, ",",
							num_test_assortments, ",",
							proc_id, ",",
							attempt, ",",
							"forest-cold-empty", ",",
							depth_limit, ",",
							forest_L1_train_error, ",",
							K, ",",
							forest_L1_test_error, ",",
							forest_approx_rate, ",",
							forest_elapsed_time)
		row_string = string(row_string, "\n");

		@spawnat 1 begin
			print(forest_outcsvhandle, row_string);
			flush(forest_outcsvhandle);
		end
	end

	## FOREST, WARM-START ##

	depth_list = [4]

	initial_left, initial_right, initial_product, initial_isLeaf = tcm_convertRankingToForest(N, orderings)

	for depth_limit in depth_list

		lambda, forest_left, forest_right, forest_product, forest_isLeaf, forest_L1_train_error, forest_elapsed_time = tcm_estimateForest_L1(N, train_assortments, choice_probs, time_limit, depth_limit,
																												initial_left, initial_right, initial_product, initial_isLeaf)
		K = length(lambda)

		forest_grand_predict, forest_grand_predict_by_class = tcm_predictForest(N, K, lambda, forest_left, forest_right, forest_product, forest_isLeaf, grand_assortments)
		forest_L1_test_error = sum( abs.( forest_grand_predict - grand_choice_probs  )  ) / 2^(N-1)

		forest_best_a = indmax( forest_grand_predict * revenues)
		forest_best_r = dot(grand_choice_probs[forest_best_a,:], revenues)
		forest_approx_rate = forest_best_r / grand_best_r

		row_string = string(single_data_name,  ",",
							num_train_assortments, ",",
							num_test_assortments, ",",
							proc_id, ",",
							attempt, ",",
							"forest-ranking-warm", ",",
							depth_limit, ",",
							forest_L1_train_error, ",",
							K, ",",
							forest_L1_test_error, ",",
							forest_approx_rate, ",",
							forest_elapsed_time)
		row_string = string(row_string, "\n");

		@spawnat 1 begin
			print(forest_outcsvhandle, row_string);
			flush(forest_outcsvhandle);
		end
	end

	helper_end_time = time()

	return helper_end_time - helper_start_time;

end

max_stars_repo_licenses: MIT Unlicense

max_issues_repo_licenses: MIT Unlicense

max_forks_repo_licenses: MIT Unlicense


text:
"""
    execute(cmd_or_cmds) -> (stdout, stderr, exitcode)

Run a Cmd/Cmds object, returning the stdout & stderr contents plus the exit code.

# Examples

```julia
out, err, code = execute(`who`)
```

If a pipeline is passed, it captures the stdout, stderr and the exitcode of the last command.

```julia
# query who are logged into the current server
out, err, code = execute(pipeline(`who`, `cut -d" " -f1`, `sort`, `uniq`))
```
"""
function execute(cmd::Cmd)
    # From https://discourse.julialang.org/t/collecting-all-output-from-shell-commands/15592
    out = Pipe()
    err = Pipe()

    process = run(pipeline(ignorestatus(cmd), stdout=out, stderr=err))
    close(out.in)
    close(err.in)
    # to avoid deadlock situation inherent in doing IO operations sequentially
    return (
        stdout = String(read(out)),
        stderr = String(read(err)),
        code = process.exitcode
    )
end

function execute(cmds::Base.OrCmds)
    out = Pipe()
    err = Pipe()

    process = run(pipeline(ignorestatus(cmds), stdout=out, stderr=err))
    close(out.in)
    close(err.in)

    # to avoid deadlock situation inherent in doing IO operations sequentially
    return (
        stdout = String(read(out)),
        stderr = String(read(err)),
        code = last(process.processes).exitcode
    )
end

"""
    on_success(f, cmd)

Apply function `f` to the stdout of shell command `cmd`.

# Examples

```julia
on_success(pipeline(`who`, `cut -d" " -f1`, `sort`, `uniq`)) do out
    filtersplit(out, "\n")
end
```
"""
function on_success(f, cmd::Union{Cmd, Base.OrCmds})
    out, err, exitcode = execute(cmd)
    if exitcode != 0
        @error err
        error("None zero exitcode $exitcode when executing `$cmd`.")
    end
    return f(out)
end


text:
using Test
using LinearAlgebra
include("../src/SSSTensor.jl")

 n = 5
 ord = 5
 tol = 1e-9
 nnz = 10
 tmppath = "tempLargestComp"

function set_up(n,ord,nnz)
    x = rand(n)
    indices = rand(1:n,nnz,ord)
    for i = 1:nnz
        indices[i,:] = indices[i,sortperm(indices[i, :])]
    end
    vals = rand(nnz)
    x, indices, vals
end

@testset "HyperGraphAlgos Tests" begin

#=
    @testset "Largest Connected Component Tests" begin

#        _, indices, vals = set_up(n,ord,nnz)
#        indices = cat(indices,(n+1)*ones(Int,ord)',1)
        indices = [1 1 2 3; 1 2 2 2; 1 3 3 3;4 4 4 4] #fully connected
        vals = [1.0,2.0,3.0,4.0]

        A = ssten.SSSTensor([(indices[i,:],vals[i]) for i in 1:nnz])
        sub_tensor,lcc_indices = get_largest_component(A::SSSTensor,filepath=tmppath)

    end
=#
    #do we need order tests?
    @testset "Connected Components Tests" begin
        indices = [1 1 2 3; 1 2 2 2; 1 3 3 3] #fully connected
        vals = [1.0,2.0,3.0]

        DictA = ssten.SSSTensor([(indices[i,:],vals[i]) for i in 1:3],3)
        COOA  = ssten.COOTen(indices,vals,3)

        assignment, sizes, _ = ssten.connected_components(DictA)
        @test assignment == ones(3)
        @test sizes[1] == 3

        COOA  = ssten.COOTen(indices,vals,3)
        assignment, sizes, orders = ssten.connected_components(COOA)
        @test assignment == ones(3)
        @test sizes[1] == 3

    end

    @testset "find_edge_incidence Tests" begin

         indices = [1 1 2; 1 2 2; 3 3 3]
         vals = [1.0,2.0,3.0]
         edges = [(indices[i,:],vals[i]) for i in 1:3]

         DictA = ssten.SSSTensor(edges,n)
         COOA  = ssten.COOTen(indices,vals,n)
         DictA_incidence = ssten.find_edge_incidence(DictA)
         COOA_incidence = ssten.find_edge_incidence(DictA)

         @test COOA_incidence == DictA_incidence
         @test COOA_incidence[1] == COOA_incidence[2]
         @test COOA_incidence[3] == Set([edges[3]])

    end

end


text:
module FESpacesTests

using Gridap
using Gridap.Algebra
using Gridap.FESpaces
using GridapDistributed
using PartitionedArrays
using Test

function main(parts)
  main(parts,SubAssembledRows())
  main(parts,FullyAssembledRows())
end

function main(parts,das)
  output = mkpath(joinpath(@__DIR__,"output"))

  domain = (0,4,0,4)
  cells = (4,4)
  model = CartesianDiscreteModel(parts,domain,cells)
  Ω = Triangulation(model)
  Γ = Boundary(model)

  u((x,y)) = x+y
  reffe = ReferenceFE(lagrangian,Float64,1)
  V = TestFESpace(model,reffe,dirichlet_tags="boundary")
  U = TrialFESpace(u,V)
  @test get_vector_type(V) <: PVector
  @test get_vector_type(U) <: PVector

  free_values = PVector(1.0,V.gids)
  fh = FEFunction(U,free_values)
  zh = zero(U)
  uh = interpolate(u,U)
  eh = u - uh

  uh_dir = interpolate_dirichlet(u,U)
  free_values = zero_free_values(U)
  dirichlet_values = get_dirichlet_dof_values(U)
  uh_dir2 = interpolate_dirichlet!(u,free_values,dirichlet_values,U)

  uh_everywhere = interpolate_everywhere(u,U)
  uh_everywhere_ = interpolate_everywhere!(u,free_values,dirichlet_values,U)
  eh2 = u - uh_everywhere

  uh_everywhere2 = interpolate_everywhere(uh_everywhere,U)
  uh_everywhere2_ = interpolate_everywhere!(uh_everywhere,free_values,dirichlet_values,U)
  eh3 = u - uh_everywhere2

  dΩ = Measure(Ω,3)
  cont  = ∫( abs2(eh) )dΩ
  cont2  = ∫( abs2(eh2) )dΩ
  cont3  = ∫( abs2(eh3) )dΩ
  @test sqrt(sum(cont)) < 1.0e-9
  @test sqrt(sum(cont2)) < 1.0e-9
  @test sqrt(sum(cont3)) < 1.0e-9

  # Assembly
  Ωass  = Triangulation(das,model)
  dΩass = Measure(Ωass,3)
  dv = get_fe_basis(V)
  du = get_trial_fe_basis(U)
  a(u,v) = ∫( ∇(v)⋅∇(u) )dΩass
  l(v) = ∫( 0*dv )dΩass
  assem = SparseMatrixAssembler(U,V,das)

  data = collect_cell_matrix_and_vector(U,V,a(du,dv),l(dv),zh)
  A1,b1 = assemble_matrix_and_vector(assem,data)
  x1 = A1\b1
  r1 = A1*x1 -b1
  uh1 = FEFunction(U,x1)
  eh1 = u - uh1
  @test sqrt(sum(∫( abs2(eh1) )dΩ)) < 1.0e-9

  map_parts(A1.values, A1.rows.partition, A1.cols.partition) do mat, rows, cols
     @test size(mat) == (num_lids(rows),num_lids(cols))
  end

  writevtk(Ω,joinpath(output,"Ω"), nsubcells=10,
    celldata=["err"=>cont[Ω]],
    cellfields=["uh"=>uh,"zh"=>zh,"eh"=>eh])

  writevtk(Γ,joinpath(output,"Γ"),cellfields=["uh"=>uh])

  A2,b2 = allocate_matrix_and_vector(assem,data)
  assemble_matrix_and_vector!(A2,b2,assem,data)
  x2 = A2\b2
  r2 = A2*x2 -b2
  uh = FEFunction(U,x2)
  eh2 = u - uh
  sqrt(sum(∫( abs2(eh2) )dΩ)) < 1.0e-9

  op = AffineFEOperator(a,l,U,V,das)
  solver = LinearFESolver(BackslashSolver())
  uh = solve(solver,op)
  eh = u - uh
  @test sqrt(sum(∫( abs2(eh) )dΩ)) < 1.0e-9

  data = collect_cell_matrix(U,V,a(du,dv))
  A3 = assemble_matrix(assem,data)
  x3 = A3\op.op.vector
  uh = FEFunction(U,x3)
  eh3 = u - uh
  sqrt(sum(∫( abs2(eh3) )dΩ)) < 1.0e-9

  A4 = allocate_matrix(assem,data)
  assemble_matrix!(A4,assem,data)
  x4 = A4\op.op.vector
  uh = FEFunction(U,x4)
  eh4 = u - uh
  sqrt(sum(∫( abs2(eh4) )dΩ)) < 1.0e-9

  dv = get_fe_basis(V)
  l=∫(1*dv)dΩass
  vecdata=collect_cell_vector(V,l)
  assem = SparseMatrixAssembler(U,V,das)
  b1=assemble_vector(assem,vecdata)
  @test abs(sum(b1)-length(b1)) < 1.0e-12

  b2=allocate_vector(assem,vecdata)
  assemble_vector!(b2,assem,vecdata)
  @test abs(sum(b2)-length(b2)) < 1.0e-12

end

end # module


text:
using CairoMakie
using naumer_Dimensionality_2022
using LinearAlgebra: tr, svd, dot, det, Diagonal, pinv, eigen
using CSV
using Random: MersenneTwister

################################################
########### Plot Parameters ####################
################################################

mrng = MersenneTwister(1234)

max_steps = 1000
state = 2*randn(mrng, 2)
state_cp = copy(state)
crlb = [4.0 0; 0 4.0]

b = -1
a = 1
τ = 0.01
γ = 0.99
actionSpace = [[sin(θ), cos(θ)] for θ in 0:.1:π]

λ = 1
psdSampleCount = [30]
trajectorySampleCount = 30
timestepSampleCount = 4
σ² = 1
max_iterations = 100
d_max = 2

T = 30

#################################################
################# Data Generation ###############
#################################################

system = HopfSystem(a, b)

optimal_sampling_trace     = zeros(max_steps)
random_sampling_trace      = zeros(max_steps)
approx_sampling_trace      = zeros(max_steps)
approx_sampling_trace_ekf  = zeros(max_steps)
optimal_sampling_trace_ekf = zeros(max_steps)

optimal_sampling_trace[1] = 8
random_sampling_trace[1]  = 8
approx_sampling_trace[1] = 8
approx_sampling_trace_ekf[1] = 8
optimal_sampling_trace_ekf[1] = 8

@info("Approximating Value Function")

values, psdSamples, stateSamples = ValueFunctionApproximation_LocalAverage_precompute(  mrng,
                                                                                system,
                                                                                τ,
                                                                                γ,
                                                                                actionSpace,
                                                                                λ,
                                                                                psdSampleCount[1],
                                                                                trajectorySampleCount,
                                                                                timestepSampleCount,
                                                                                σ²,
                                                                                max_iterations,
                                                                                d_max)

@info("Solving Optimal Sampling Dynamic Programming Problem")

NN_Policy_tuple = (system, σ², τ, values, psdSamples, stateSamples, actionSpace, d_max)

for i in 2:max_steps
    global state
    global crlb
    global NN_Policy_tuple

    # Choose Action and Update System
    action, index, state, crlb = LocalAverage_OptimalPolicy(    state,
                                                                crlb,
                                                                NN_Policy_tuple...)

    # Record cost
    optimal_sampling_trace[i] = tr(crlb)
end

@info("Solving Optimal Sampling Dynamic Programming Problem - EKF")

NN_Policy_tuple = (system, σ², τ, values, psdSamples, stateSamples, actionSpace, d_max)
beliefState = copy(state_cp) + 2*randn(2)
beliefCRLB = [4.0 0; 0 4.0]
crlb = [4.0 0; 0 4.0]

for i in 2:max_steps
    global state
    global crlb
    global beliefState
    global beliefCRLB

    # Choose Action
    action, ~, ~, ~ = LocalAverage_OptimalPolicy(   state,
                                                    crlb,
                                                    NN_Policy_tuple...)

    # Update True System
    jacobian = flowJacobian(state, τ, system)
    crlb = updateCRLB(crlb, action, jacobian, σ²)
    state = flow(state, τ, system)

    # Update Belief System
    observation = dot(action, state) + sqrt(σ²) * randn(mrng)
    jacobian = flowJacobian(beliefState, τ, system)
    beliefCRLB = updateCRLB(beliefCRLB, action, jacobian, σ²)
    beliefState = flow(beliefState, τ, system)
    beliefState = stateupdate_EKF(beliefState, beliefCRLB, observation, action, σ²)

    # Record cost
    optimal_sampling_trace_ekf[i] = tr(crlb)
end

@info("Solving Optimal Sampling 1D Approximation Problem")

crlb = [4.0 0; 0 4.0]
state = copy(state_cp)

for i in 2:max_steps
    global state
    global crlb
    global T
    global σ²
    global system

    # Decide action
    ~, action = optimalaction_1DApprox(actionSpace, system, state, crlb, σ², T)

    # Update System
    jacobian = flowJacobian(state, τ, system)
    crlb = updateCRLB(crlb, action, jacobian, σ²)
    state = flow(state, τ, system)

    # Record Cost
    approx_sampling_trace[i] = tr(crlb)
end

@info("Solving Optimal Sampling 1D Approximation Problem - EKF")

crlb = [4.0 0; 0 4.0]
state = copy(state_cp)
beliefState = state + 2*randn(2)
beliefCRLB = [4.0 0; 0 4.0]

for i in 2:max_steps
    global state
    global crlb
    global beliefState
    global beliefCRLB

    # Decide action
    ~, action = optimalaction_1DApprox(actionSpace, system, beliefState, beliefCRLB, σ², T)

    # Update True System
    jacobian = flowJacobian(state, τ, system)
    crlb = updateCRLB(crlb, action, jacobian, σ²)
    state = flow(state, τ, system)

    # Update Belief System
    observation = dot(action, state) + sqrt(σ²) * randn(mrng)
    jacobian = flowJacobian(beliefState, τ, system)
    beliefCRLB = updateCRLB(beliefCRLB, action, jacobian, σ²)
    beliefState = flow(beliefState, τ, system)
    beliefState = stateupdate_EKF(beliefState, beliefCRLB, observation, action, σ²)

    # Record Cost
    approx_sampling_trace_ekf[i] = tr(crlb)
end

@info("Generating Random Sampling Trajectory")

crlb = [4.0 0; 0 4.0]
state = copy(state_cp)

for i in 2:max_steps
    global state
    global crlb

    # Choose Action
    action = actionSpace[rand(mrng, 1:length(actionSpace))]

    # Update System
    jacobian = flowJacobian(state, τ, system)
    crlb = updateCRLB(crlb, action, jacobian, σ²)
    state = flow(state, τ, system)

    # Record Cost
    random_sampling_trace[i] = tr(crlb)
end

#################################################
################ Save Data CSV ##################
#################################################

CSV.write("HopfSampling.csv", ( random_sampling_trace = random_sampling_trace,
                                optimal_sampling_trace = optimal_sampling_trace,
                                optimal_sampling_trace_ekf = optimal_sampling_trace_ekf,
                                approx_sampling_trace = approx_sampling_trace,
                                approx_sampling_trace_ekf = approx_sampling_trace_ekf,
                                samples = Array(1:max_steps)))


text:
using DataFrames, CSV, Statistics, Dates, Plots, Measures; pyplot()

data = CSV.read("data/temperatures.csv", DataFrame)
brisbane = data.Brisbane
goldcoast = data.GoldCoast

diff = brisbane - goldcoast
dates = [Date(
            Year(data.Year[i]),
            Month(data.Month[i]),
            Day(data.Day[i])
        ) for i in 1:nrow(data)]

fortnightRange = 250:263
brisFortnight = brisbane[fortnightRange]
goldFortnight = goldcoast[fortnightRange]

default(xlabel="Time", ylabel="Temperature")
default(label=["Brisbane" "Gold Coast"])

p1 = plot(dates, [brisbane goldcoast],
        c=[:blue :red])
p2 = plot(dates[fortnightRange], [brisFortnight goldFortnight],
        c=[:blue :red], m=(:dot, 5, Plots.stroke(1)))
p3 = plot(dates, diff,
        c=:black, ylabel="Temperature Difference",legend=false)
p4 = histogram(diff, bins=-4:0.5:6,
        ylims=(0,140), legend = false,
        xlabel="Temperature Difference", ylabel="Frequency")
plot(p1,p2,p3,p4, size = (800,500), margin = 5mm)

max_stars_repo_name: newptcai/StatsWithJuliaBook

max_issues_repo_name: newptcai/StatsWithJuliaBook

max_forks_repo_name: newptcai/StatsWithJuliaBook


text:
module TestThis

# package code goes here
println("Sherlock! We are loading the module now...")
using PyPlot

function greet(name)
	println("Hello, $name")
end

function addme(a,b)
	return a + b
end

function plotStuff(savehere::String)
	a = [1,2,3,4]
	b = [5,4,5,3]
	plot(a,b)
	savefig(savehere*".png")
end

end # module


text:
using DynamicGrids, Test
using DynamicGrids: isshowable, allocateframes!, frameindex, storeframe!, SimData

init = [10.0 11.0;
        0.0   5.0]

output = ArrayOutput(init)
ruleset = Ruleset(Life())

@test frameindex(output, 5) == 5
@test isshowable(output, 5) == false

# Test pushing new frames to an output
update = [8.0 15.0;
          2.0  9.0]
@test length(output) == 1
push!(output, update)
@test length(output) == 2
@test output[2] == update

# Test allocateing additional frames to an output
allocateframes!(output, init, 3:5)
@test length(output) == 5
@test firstindex(output) == 1
@test lastindex(output) == 5
@test size(output) == (5,)

# Test storing a new frame
@test output[3] != update
data = SimData(ruleset, update, 1)
storeframe!(output, data, 3)
@test output[3] == update

# Test creting a new output from an existing output
output2 = ArrayOutput(output)
@test length(output2) == 5
@test output2[3] == update